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Abstract 



The present work investigates the numerical evolution of linearized oscillations of non-rotating, spheri- 
cally symmetric neutron stars within the framework of general relativity. 

We derive the appropriate equations using the (3+l)-formalism. We first focus on the evolution of 
radial oscillations, which do not emit gravitational waves. The associated system of equations being 
quite simple, we demonstrate how to handle a numerical instability that also occurs in the non-radial 
case, when the stellar model is constructed based on a realistic equation of state. We devise a coordinate 
transformation that not only removes this instability but also provides much more accurate results. For 
comparison reasons, we compute the eigenfrequencies of the radial modes with an eigenvalue code and 
thereby confirm the results of Vath & Chanmugam, which differ from previous calculations that were 
performed by Glass & Lindblom. 

The main part deals with the evolution of non-radial oscillations (/ > 2) of neutron stars. Here, we 
compare different formulations of the equations and discuss how they have to be numerically dealt with 
in order to avoid instabilities at the origin. We present results for various polytropic stellar models and 
different initial data. They show that the quasi-normal modes of the star, such as the /-, the p-, and the w- 
modes, can, indeed, be excited by suitable initial data. However, the excitation strength of the it;-modes 
strongly depends on the chosen initial data. For some initial data the occurrence of the u)-modes can be 
totally suppressed. 

For ultra-compact models we find the interesting feature that the first ring-down phase of the wave 
signal cannot be associated with any of the known quasi-normal modes that belong to the star itself; the 
frequency and damping time rather correspond to the first quasi-normal mode of an equal mass black 
hole. 

When switching to realistic equations of state, we find that we face the same numerical problems as in 
the radial case. Here, too, we can get rid of them by means of the same coordinate transformation inside 
the star, but things are more comphcated because the fluid equation is coupled to the metric equations, 
which propagate the gravitational waves. For those equations the transformation is not defined, therefore 
we have to interpolate between the different grids. For polytropic equations of state this can give rise to a 
quite strong violation of the Hamiltonian constraint, which can spoil the resulting wave signal. However, 
for realistic equations of state, and it is only for those that the transformation is necessary, this does not 
happen and the equations can be integrated in a stable way and also provide quite accurate results. 

In the last part of this thesis we consider a physical mechanism for exciting oscillations of neutron 
stars. We use the time dependent gravitational field of a small point mass /i that orbits the neutron star 
to induce stellar oscillations. Hereby, we assume /i to be much smaller than the mass M of the neutron 
star. In this so-called particle limit, the gravitational field of the moving particle is considered to be a 
perturbation of the background field of the neutron star. With this particle we have a physical means 
which removes the arbitrariness in choosing the initial data. 

However, even with the presence of the particle, there is still too much freedom in constructing the 
initial data, which is due to the fact that in addition to the field of the particle we always can superpose 
any arbitrary amount of gravitational waves. Our task is to minimize this additional gravitational- wave 
content and find initial data which correspond to the pure gravitational field of the particle. By looking at 
the flat space case we can construct analytic initial data that satisfy the above requirements. Those then 
will serve as a good approximation for the "real" initial data in the curved spacetime of the neutron star. 
By sampling various orbital parameters of the particle we show that in general the particle is not able to 
excite any u;-modes. It is only for speeds very close to the speed of hght that the u;-mode is a significant 
part of the wave signal. This result indicates that it is rather improbable that any physical mechanism 
which can be simulated by an orbiting particle can excite the u;-modes of neutron stars in a significant 
manner. 
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Chapter 1 
Introduction 



The investigation of stellar oscillations has a quite long history, and is a whole astrophysical 
branch on its own. The evaluation of the oscillation modes of stellar objects has helped to 
reveal a lot of information about their interior structure. With the use of asteroseismological 
methods we now have a quite detailed understanding of the physics of a whole range of stellar 
objects, be it our own sun or a white dwarf somewhere in our galaxy. 

Most stellar objects can be adequately described by Newtonian theory, however, for very 
compact objects such as neutron stars, the effects of general relativity cannot be longer ne- 
glected. In fact, they are very important, for now, stellar oscillations will be associated with the 
emission of gravitational waves, which can carry important information about the physics inside 
those compact objects. These waves, when detected, will open a totally new observational win- 
dow for astrophysicists. Gravitational waves are not affected by events at the surface of the star, 
which is what happens to electromagnetic radiation and they also remain practically unaffected 
by any kind of matter while travelling through space. Thus, they carry "clean" information of 
the physical properties of the neutron stars. 

It is well known that lowest gravitational-wave multipole is the quadrupole radiation, hence 
the radial and dipole oscillations do not emit gravitational waves. From the gravitational-wave 
astronomical point the latter are therefore quite uninteresting, however, they could make them- 
selves visible through tiny undulations in the electromagnetic radiation signal of the neutron 
stars. 

The study of the non-radial oscillations modes of non-rotating neutron stars within the 
framework of general relativity was initiated in 1967 by a series of papers by Thorne and sev- 



eral coauthors [jrg, [T^, |2^, |2T], |23]. In the following decades a lot of authors made endeavors to 
study the relativistic oscillation spectrum of neutron stars, which turned out to be particularly 
rich. 

From the general Newtonian theory of non-radial linear oscillations it follows that the oscil- 
lations of a spherical stellar object can be divided into two classes according to their transfor- 
mation behavior under space reflection. One class consists of even parity or polar modes; the 
odd parity or axial modes belong to the other class. 

If the neutron star is modeled with a perfect fluid then from Newtonian theory it follows that 
any non-radial oscillations such as the /undamental mode [BO], the pressure modes, and the 
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(gravitational modes must belong to the polar class. This means that for Newtonian perfect 
fluid stellar models there are no axial modes at all. It is only with the inclusion of a solid crust 
that there can exist axial torsional modes. 

In the first 20 years after the fundamental paper of Thome, it was common belief that the 
oscillation spectra of relativistic neutron stars would not differ much from Newtonian stars. The 
only new effect of general relativity was thought to be the radiation-damping of the oscillation 
modes due to the emission of gravitational waves. Therefore the main focus was to compute the 
frequencies and damping times of those polar modes which also exist in the Newtonian theory, 
but now within the framework of general relativity. The only axial modes under consideration 
were the t-modes of neutron stars with a solid crust [52l |3l, HO]. 



It was only in 1988 that another family of polar modes was found by Kojima [|38|], whose 
existence has already been anticipated two years earlier in a toy model by Kokkotas & Schutz 
[^. Those new modes, which were termed wave-modes by Kokkotas & Schutz [|45|], are of 
purely relativistic origin and do not have a Newtonian counterpart. They are predominantly 
metric oscillations and couple only weakly to the motion of the matter. 

Around the same time Chandrasekhar & Ferrari showed that for ultra-relativistic stellar 
models, there exists also a family of axial perturbations, which do not couple to the matter at all. 
This is because the spacetime curvature inside the star can be so strong that it can trap impinging 
gravitational waves. Those "trapped" modes are quite long-lived since they correspond to quasi- 
bound states inside the gravitational potential of the neutron star, which only slowly leak out. 

A little later, it was found by Leins et al. [ p8| ] that the polar w-modes are split into two 
branches, one with only a few strongly damped modes (called wn- or interface modes), and the 
other one with an infinite number of modes with increasing frequency and damping times. 

In 1994 Kokkotas [ pO[ ] showed that for less relativistic stars, which do not possess any 
trapped modes, there also exists a family of axial modes, which are quite similar to the polar w- 
modes. Finally, in 1996 Andersson et al. [ pT] , p2| ] established that all three kinds of gravitational- 
wave modes (the iw-modes, the w//-modes and the trapped modes) exist for both the polar and 
the axial cases, the trapped modes, however, only for ultra-relativistic star. 

This rich set of various pulsation modes having been found, it then was natural to ask 
whether they would indeed be excited in a real pulsation scenario. Allen et al. [ |5^ took a 
first step in answering this question by numerically integrating the evolution equations for sev- 
eral different sets of initial data. By looking at the emitted gravitational radiation, they indeed 
found that both the fluid /- and p-modes and the gravitational w-modes can be excited. How- 
ever, the excitation strength depends strongly on the choice of initial data, which, of course, 
were constructed ad hoc and had no astrophysical meaning whatsoever. 

Andersson & Kokkotas [ ]55| ] have shown that, once the parameters of the /- and the first 
w-mode are known, it is in principle possible to obtain the mass and the radius of a neutron 
star quite accurately. Those quantities could then be used to determine the equation of state 
for the neutron star matter. Because of the still quite poor understanding of the physics at the 
subnuclear density level, there exists a plethora of different equations of state in the literature, 
all based on different physical models of the nuclear interaction. But once the mass and radius 
of a neutron star are known with a certain accuracy, it is possible to rule out those equations of 
state that predict a different radius for a given mass, or vice versa. And this in turn would hint 
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at the physics that governs the behavior of the neutron star matter in the center of the star. 

It is therefore of particular interest whether or not the w-modes will be excited in a real 
astrophysical scenario such as the collapse of a progenitor star to a neutron star. Whereas the 
/- and some fluid p-modes will always be present in an wave signal, it is not clear if this is also 
the case for the u'-modes. The numerical experiments show that it is possible to construct initial 
data which do not excite any w-modes at all. 

To satisfactorily answer the question what a real signal would look like, one would have 
to follow the appropriate scenario (e.g. the collapse of a burned out star core) in full detail 
from the point on where relativistic effects become important. Unfortunately, this is, for the 
time being, still impossible. This is mainly because one would have to use the full set of the 
nonlinear Einstein equations together with the whole machinery of nuclear physics to simulate 
a core collapse. However, when the newborn neutron star has almost settled, the study of stellar 
pulsations can be done within the linearized theory, which makes it tractable. But even there 
one usually restricts oneself to very simple neutron star models. In our case, we will neglect the 
spin of the star, which makes the equations much simpler, but also prevents us from studying 
the effects of rotation, which might be very important. 

There are no unstable non-radial modes for non-rotating stars as long as the Schwarzschild 
discriminant 

S(r) = ^-(^] ^ 
dr \de J ^ dr 

is positive [p7|, [7T|]. For negative S the oscillations become unstable with respect to convection. 
In this thesis we always use barotropic equations of state, for which it is 5 = 0, therefore all 
non-radial modes will be stable. 

However, things can change dramatically when rotation is included. There will be always 
some modes which will undergo the CFS -instability |59| , ^], and recently a whole set of 
unstable modes (r-modes) has been discovered. They might be responsible for slowing down 



rapidly rotating neutron stars [p6|]. For a recent review of the stability properties of rotating 
neutron stars, see [[7l|]. Other recent reviews on oscillations of neutron stars and black holes can 
be found in [|5|,|6^]. 

The thesis is organized as follows: 

In chapter 2 we will derive the general form of the perturbation equations using the (3-1-1)- 
split. 

In chapter 3 we will specialize on the radial case, where we investigate both the evolution 
and the eigenvalue problem. We discuss the occurrence of a numerical instability, which occurs 
both in the radial and non-radial cases when the neutron star is modeled with a realistic equation 
of state. We also present a remedy which not only removes the instability but also automatically 
increases the accuracy of the evolution as compared to the old set of equations. 

In chapter 4 we focus on the non-radial oscillations. We present various forms of the equa- 
tions and discuss which of them are the most suitable set for the numerical evolution. We show 
how to treat the boundary conditions in order to obtain stable evolutions for any value of /. 
Here, too, we demonstrate how to get rid of the same numerical instability that occurs with the 
use of realistic equations of state. 
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In chapter 5 we examine the excitation of neutron star oscillations by the scattering of a 
small point mass fj, in the particle limit, where we assume fx to be much smaller than the mass 
M of the neutron star. The gravitational field of the particle, which follows a geodesic path in 
the background metric of the neutron star, can then be treated as a perturbation of the back- 
ground. We sample a variety of orbital parameters to see whether or not the particle is able to 
excite the w-modes of the neutron star. 



Notations and conventions 



• We use Einstein's sum convention. Greek indices run from to 3, Latin indices from 1 to 
3. The time component is the 0-component. 

• The signature of the metric is (— , +, +, +), i.e. time-like 4-vectors have negative norm. 

• We use geometric units, i.e. we set c = G = 1. 

• Covariant derivatives with respect to are denoted by and partial derivatives ^ are 
sometimes abbreviated by 9^. 

• Derivatives with respect to the radial coordinate r are often denoted by a prime, and time 
derivatives by an over-dot. 

• The perturbation variables should be properly denoted by Qim{'t,i^)'> however, we will 
often omit the indices / and m. 



Chapter 2 

Linearizing Einstein's equations 



Within perturbation theory we first construct a stationary non-pulsating and non-rotating neu- 
tron star model by solving Einstein's equations, which in this case reduce to a set of coupled 
ordinary differential equations. 

Having found this (numerically) exact background solution g^^,^, we describe the pulsations 
by "small" deviations from the original metric. Thus, the metric of the perturbed spacetime 
will be written as g^^, — g^,^ + h^i,, where h^j,^, is considered to be a small perturbation of the 
background metric g^i,. 

To find the relevant equations for h^y, we plug the perturbed metric g^i, into Einstein's 
equations and neglect all terms that are quadratic or of a higher power in the h^^. We will thus 
obtain a linear system of partial differential equations for the h^^, which is much easier to solve 
than the full nonlinear set. Still, this set of equations depends on all four coordinates, which, 
even with the presently available computer power, would be too time-consuming to be solved 
in its fuUfledged form. 

However, because of the spherical symmetry of the background metric we can get rid of the 
angular dependence by expanding the perturbation equations into spherical tensor harmonics. 
It is thus possible to reduce the equations into a (l+l)-dimensional evolution system, which can 
be numerically solved on present day PCs. 

2.1 The unperturbed stellar model 

The simplest model for a neutron star is a non-rotating zero temperature perfect fluid sphere, 
whose static spherically symmetric geometry is given by a line element of the form 

ds"" = -e^^df + e'^^dr^ + r^dO^ + sin^ ed(j)^) , (2.1) 

where the two functions v and A only depend on the radial coordinate r and have to be deter- 
mined by the field equations. The appropriate energy-momentum tensor is given by 

T^,^ = (e + p) u^u,, + pg^,, , (2.2) 

where p is the pressure, e the energy density, and the covariant 4-velocity of the fluid. In the 
rest frame of the fluid, which is static, the only non-vanishing component is the time component 
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uo = —e", and Einstein's equations Gfj,u = SnT^iy and the conservation equations D^T^" = 

yield the following three structure equations for the four unknown A, p, and e: 

1 - 

A' = + 47rre2^e (2.3a) 

2r 

u' = + 4:nre^^p (2.3b) 

2r 

p' = — z/'(p + e). (2.3c) 

To fully determine this system of equations, an equation of state 

p = p{s,n) (2.4) 

e = e(s,n) (2.5) 



must be supplemented. We will refer to the equations ( [2 .SI ) as Tolman-Oppenheimer-Volkov or 
TOV equations, even if the original TOV equations are written in a slightly different form. 

Throughout this work we always assume the neutron star to have zero temperature, which 
is quite reasonable since the pressure inside the neutron star is mainly maintained by a Fermi 
gas of degenerate neutrons. Hence, the specific entropy s can be set to zero, too, and we can 
eliminate the baryon density n and obtain a barotropic equation of state, where the pressure is a 
function of the energy density alone: 

p = p{e) . (2.6) 

The simplest and therefore quite often used form is given by a so-called polytropic equation of 
state 

P = «:e^ (2.7) 

where k and V are the polytropic constant and polytropic index, respectively. More realistic 
equations of state have to include the microphysics that dictates the interplay between p and e 
on the nuclear and subnuclear levels for the neutron star matter. However, the physics under 
the extreme conditions of the high pressures that prevail in the center of neutron stars is not yet 
fully understood. In the literature there are therefore quite a few different realistic equations 
of state, which were calculated based on various different microscopic models of (sub)nuclear 
interactions. 

It is clear that the zero temperature assumption neglects all kinds of thermal and viscous ef- 
fects which can affect the stellar oscillations. On the one hand it will suppress the existence of 
the whole family of ^f-modes, and on the other hand we ignore the effects of viscosity and inter- 



nal friction, which would normally damp out any oscillation. However, the damping times [ pTQ 
are such that their neglect will have no effect on the numerical evolutions we are investigating. 
If we introduce an additional function m, which is related to A by 



-2A ^ 2m 



e-'^ = 1 , (2.8) 



2.2. The perturbation equations 
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we obtain from ( [2.3a| ) 

m' = ATTr^e . (2.9) 
Integration immediately leads to 

m(r) = / Arrf'^edf , (2.10) 
Jo 

which makes it clear that m(r) represents the total gravitational mass enclosed inside the radius 
r. 

To obtain the stellar model, we have to integrate the TOV equations ( |2.3| ) together with ( |2.6D 
from the center up to the point where the pressure p vanishes. This then defines the surface R 
of the star. 

From Birkhoff theorem it follows that the exterior vacuum region of the star is described by 
the Schwarzschild metric with the mass parameter M = m{R). 



2.2 The perturbation equations 



A somewhat more elaborate derivation of the relevant equations can be found in my Diploma 



thesis PSBQ. Here, we will only briefly discuss the necessary steps and not go into very great 
details. Our starting point is the ADM-formulation [ [73| ] or (3+l)-decomposition of the Einstein 
equations. In this formulation the 10 field equations are first split into 6 dynamical equations 
and 4 constraint equations. The dynamical equations, which have second order time derivatives, 
are then cast into a set of 12 evolution equations, which are first order in time. Those are the 
6 evolution equations for the 3-metric 7ij of a space-like 3-dimensional hypersurface S and 
another 6 equations for the time development of its extrinsic curvature Kif 



dtKij 



a 



(2.11) 



(2.12) 



Herein, a denotes the lapse function and is the shift vector. The remaining 4 constraint 
equations, which have to be satisfied by any physically acceptable initial data, are given by 



R - KijK'^ + 
DiK - DjK\ 



IQnp 
8nji . 



(2.13) 
(2.14) 



Here, p and ji are the energy density and momentum density measured by a momentarily sta- 
tionary observer, whose 4-velocity coincides with the time-like normal vector of the space- 
like hypersurface: 



(2.15) 
(2.16) 
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The scalar equation ([2.13P is usually called the Hamiltonian constraint and the vector equation 



( [2.14D is called momentum constraint. Once satisfied on the initial hypersurface, they will be 



automatically preserved throughout the evolution by virtue of the Bianchi-identities. 

As we have already mentioned above, the oscillations of the neutron star will be described 
within the framework of perturbation theory, i.e. we will treat them as small perturbations 
around the fixed background which is given by the unperturbed stellar model. Hence, we 
write the perturbed metric g^y as a sum of the static background g^^ and the time-dependent 
perturbations h^^: 

Qi^v = 9tiu + h^,y . (2.17) 



Since we are using the (3+l)-decomposition, we have to write down the perturbed metric ( |2.17[ ) 



in terms of lapse function, shift vector, and 3-metric. Furthermore, we have to deal with the 
extrinsic curvature as a dynamical variable. Their respective perturbations will be denoted by 



a, P\ hij, and kij. For the background metric (|2.1D, shift and extrinsic curvature vanish and the 
unperturbed lapse function A is given by 



A = = e^ (2.18) 

In addition to the metric perturbations, we have to describe the perturbations of the energy- 
momentum tensor T^^. For a perfect fluid, the only quantities that can be perturbed are energy 
density e, pressure p, and 4- velocity u^j, whose perturbations will be denoted by de, 6p, and 
respectively. For a barotropic equation of state, however, 5e and 5p are not independent but are 
rather related through 

Sp = ^Se, (2.19) 
ae 

where ^ =: is the square of the sound speed Cg inside the fluid. 



Due to the simple form of (|2.lD, the evolution equations for the 3-metric and the extrinsic 



curvature perturbations do not become as messy as they usually do in perturbation theory: 

dthij = + ikidjf3'' + 7fc,9./5'= - 2e^% (2.20) 

(2.21) 



+ [5Ri, + An {{p - e)h,, + 5eiC^^ - 1)7.,)] • 

Herein, 'jij denotes the spatial part of (P!T|), its inverse is 7*^, and the perturbed Christoffel 
symbols are defined as 

= ^7"" {d^hn^j + d^h^i - dmk, - 2r\^hi^) . (2.22) 
The perturbed Ricci tensor is given by 

6R,j = D,6T% - Dj6T\, 

= dk6T% - d,6T',, + T\^6T'i, + T),ST\^ - T\,6T% - T)^6T\, . (2.23) 



2.2. The perturbation equations 
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To first order in the perturbations, the constraints ( |2.13| ) and ( [2.14[ ) read 

Y^dRij-h'^ij = lQiT6e (2.24) 
yfc ^0,k,k - d,hk - V\^k,i + r^.A;,,) = -87r(p + e)5ui . (2.25) 

It is interesting to note that it is possible to eliminate 5e in ( [2.21D by virtue of the Hamiltonian 



constraint ( [2.24| ) and to obtain thus a consistent system of evolution equations for the metric and 
extrinsic curvature alone. The constraints then can serve to compute the matter perturbations 
5e and 5ui. In this case we would not have to use the equations of motion for the matter 
perturbations, which follow from the conservation of energy-momentum D^T^^^ . 

Our next step consists of expanding the equations in spherical tensor harmonics. Any 
symmetric tensor A^^ can be expressed in terms of a set of ten spherical tensor harmonics 
{[yLVui.^^ 0)}a=i,...,io, namely 

oo I 10 

A,At, r,e,<P) = ^Ait, r) [ytUiO, 0) • (2.26) 

1=0 m=-l A=l 

Due to the spherical symmetry of the background we thereby can eliminate the angular de- 
pendence and obtain equations for the coefficients ci^(t, r). The field equations thus will be 
reduced to partial differential equations in t and r. 

There are various ways to define those tensor harmonics 3^/^^, and an exhaustive overview 
can be found in The set first given by Regge and Wheeler [jT7|] is widely used throughout 
the literature and we will follow this tradition, even if this set has the disadvantage of not being 
orthonormal. However, this is of no account in the following proceeding. It is only with the 
inclusion of an orbiting particle that this fact causes some minor inconveniences. 

The Regge- Wheeler harmonics can be divided into two subsets that behave in different ways 
under parity transformation. Under space reflection the polar or even parity harmonics change 
sign according to (—1)', whereas the axial or odd parity harmonics transform like (—1)'+^. 

We now expand the metric as follows: 
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Im 

+ vny'LU + vt'iyLU + vfiyLU + vtiyLU (2-27) 



+ Tt'iyLu + TtiyLu + T^yi 



The notation has been chosen such that the coefficients Si represent the scalar parts of h^^, 
namely a,l3r, and h„, whereas the Vi stand for the vector components Pe^P^.h^e and /i^^. 
Lastly, the % represent the tensorial components hee, Iiq^ and h^^f,. Note that this expansion 
includes both polar and axial harmonics. Similarly, the extrinsic curvature tensor kij will be 
expanded as 



= Ki"^[yLh + K^2"'[yLh + KinyJ 



+ Kt[yLh + Kt[yLh + Ktiy'' 



^n^J^. (2.28) 



lm\ij 
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Last not least, we need the matter variables 

de = p'^^Yi^ (2.29) 

= «f [3^L].i + u'nyLu + n'nyJmU ■ (2.30) 



Here, Yim in the expansion ( |2.29D for 5e is just the ordinary spherical harmonic. Using this 
expansion, we will obtain 12 evolution equations for the coefficients of the 3-metric hij and the 
extrinsic curvature kij. However, it is clear that this set cannot be used for the evolution, for we 
have not specified any gauge yet. Picking a specific gauge means prescribing the coefficients 
of lapse a and shift I3i for which we do not have any evolution equations. Of course, there are 
various ways to choose lapse and shift, and for the numerical evolution of the full nonlinear 
Einstein equations it is crucial to make a good choice because inadequate lapse or shift may 
cause the code to crash at very early times due to the blow up of some variables. 

However, we are dealing with the linearized field equations and we do not have to worry too 
much about some of the nasty things that might happen in the nonlinear case. We will choose 
lapse and shift in a way as to simplify the equations as much as possible. 

We could now proceed and write down the general form of the perturbation equations, but 
we will not do so because on the one hand the equations are quite lengthy and on the other hand 
they are not very useful as such. Besides, they can be found in Appendix C of my diploma 



thesis PSBQ. However, it seems to be appropriate to make some comments on the equations. 

First of all, the equations are independent of the spherical harmonic index m because of the 
spherical symmetry of the background. Secondly, they decouple with respect to their behavior 
under parity transformation. Hence we obtain two sets of equations, one describing polar pertur- 
bations, the other axial ones. The latter cannot generate any density or pressure changes in the 
neutron star for those quantities are scalars and therefore have even parity. The equations can be 
further subdivided according to their value of /. The radial and dipole modes are characterized 
by / = and / = 1, respectively. Those oscillations do not give rise to gravitational radiation, 
hence the equations are only meaningful in the interior of the star itself. Of course, they are also 
valid in the exterior, but we can make them identically vanish by means of an appropriate gauge 
transformation. Furthermore, for / = we do not have any angular dependence at all, hence all 
the metric and extrinsic curvature coefficients that are proportional to a derivative of Yim vanish 
identically. We will consider this case in the next chapter. The dipole modes, which have been 
studied in [ pS] , |T], will not be treated in this thesis. 

For / > 2 it is not possible any more to make the exterior equations vanish, hence those 
equations describe real gravitational waves, which propagate through the spacetime and carry 
information about the oscillations of the neutron star. This case is the main point of investigation 
of this thesis and is presented in chapter 4. 



Chapter 3 

Radial oscillations of neutron stars 



As they are the simplest oscillation modes of neutron stars, radial modes have been the first 
under investigation [ |T^ ] . More important, they can give information about the stability of the 
stellar model under consideration [1^]. Since they do not couple to gravitational waves, 
the appropriate equations are quite simple and it is relatively easy to numerically solve the 
eigenvalue problem that leads to the discrete set of oscillation frequencies of a neutron star. 
In the absence of any dissipative processes the oscillation spectrum of a stable stellar model 
forms a complete set; it is therefore possible to describe any arbitrary periodic radial motion of 
a neutron star as a superposition of its various eigenmodes. It hence seems quite superfluous to 
explicitly solve the time dependent equations, for we do not expect to gain new physical insight. 
This might be true indeed. There are, however, quite unexpected numerical problems that are 
associated with the evolution of the radial oscillations of realistic neutron star models. 

Those problems are quite generic, for they also occur for the non-radial oscillations and 
probably also in the case of rotating neutron stars, and once we have them under control in the 
radial case it should be straightforward to confer the appropriate numerical treatment to other 
cases. 

The above mentioned numerical problems are instabilities that occur when the neutron star 
model is constructed with a realistic equation of state. However, they do not appear, if one 
uses polytropic equations of state, which is often done for the sake of simplicity. As it turned 
out, this instability is not a general instability of the numerical scheme, for it is dependent 
on the resolution and can be made to disappear if the resolution exceeds a certain threshold. 
However, this threshold strongly depends on the numerical scheme that is used to evolve the 
equations, and can be so high that it prevents any evolutions within a reasonable time limit. For 
other discretizations, the threshold can be relatively low and does not represent a real obstacle 
to obtaining numerical results. However, this is mainly because the radial case is a (l+l)-d 
problem and one only has to consider the stellar interior, since the exterior spacetime remains 
totally unaffected. For non-radial oscillations with / > 2 this is not true any more, for here 
the oscillations will generate gravitational waves, which propagate towards infinity. This means 
that we have to include the exterior domain in our numerical evolution as well, which will result 
in a much bigger computational expenditure. Here, the required minimum resolution can extend 
the computation time of a single run to quite large values. Still, even the non-radial case is a 
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(l+l)-d problem, and if we include rotation we will face a (2+l)-d problem, which can become 
totally intractable if a high resolution is required to obtain a stable evolution. 

As we shall see, the cause of the instability is confined to a small region close to the surface 
of the star, and it is only here that the high resolution is needed in order to obtain stability. For 
the more time-consuming cases, we then could only locally refine the grid in this region instead 
of using high resolution for the whole domain. However, for each stellar model, we would 
have to localize the troublesome region and find the required resolution, which is more or less 
a matter of trial and error. 

Therefore, we seek for a better way to solve this problem and we find it in reformulating the 
equations in such a way that the instability totally disappears, and the equations can be evolved 
for any stellar model and for any desired resolution in a stable way. 



3.1 Derivation of the evolution equations 

The first one to write down the equations was Chandrasekhar |T^, later on various authors 
rewrote them in many different ways [ pT| , [T3|, |1^. Since the (3+l)-formalism is particularly 
suitable for the numerical evolution we will rederive the equations using the framework of the 
(3+ 1 )-decomposition. 

From the last chapter we know that radial perturbations are described by / = 0. Since 
the appropriate harmonic is just a number Yqo = l/v^47r there is no angular dependence at all 
and any derivative of Iqo with respect to ^ or vanishes. In this case we can absorb Fqo in 
the perturbation variables, which are then functions of t and r only, and the expansion of the 
perturbations reads 

a = e''Si{t,r) (3.1a) 

Pr = re^^S2{t,r) (3.1b) 

re2^53(t,r) 
hij = I r'^T{t,r) | . (3.1c) 

sin^e T{t,r) 

Similarly, we have for the extrinsic curvature 

e^^Ki{t,r) 
\r^K2 


This particular decomposition has been chosen in order to obtain a more convenient set of 
equations. The matter perturbations are characterized by the perturbation of the energy density 
5e and the (covariant) radial component of the 4-velocity Uj., which are expanded as 




= -e-" \r^K2{t,r) | . (3.2) 



5e = p{t,r) 
Sur = —e^uit^r) 



(3.3) 
(3.4) 



3.1. Derivation of the evolution equations 



15 



We then obtain the following four evolution equations for the metric perturbations and the 
perturbations of the extrinsic curvature: 



dt 
dT 

'dt 
dKi 



dr 

= 2S2 + K2 

- Qri.' + 1 

+ 47re^^ (1 - C',) p 



dT 



+ r-^'hz: + (2-'-A') 



dr 



or 



,2X 



3 , e 

-u' H 

2 r 



- 2 



dr 
S3 



(3.5a) 
(3.5b) 

(3.5c) 



2u-2X 



,2A 



dK2 

dt 



+ STre^'^ (1 - CI) p . 
In addition we have the Hamiltonian constraint 



r J dr r dr dr 



+ 2— T + 2A' - 2i/' - - S: 



(3.5d) 



,2A 



+ A'- - -T^ 



dr 



^^^+2(--A')53 



9r 



and the momentum constraint 

dK^ 
dr 



Stt (p + e) e^^'u 



- K2 + -K,. 
r I r 



(3.6) 



(3.7) 



Note, that in the above equations we have not yet specified any gauge. To perform numerical 
evolutions, we have to fix the gauge, that is we need some additional prescriptions for lapse Si 
and shift 5*2. Let us choose vanishing shift 



-52 = 0, 



(3.8) 



but, for the moment, let us keep the lapse still unspecified. Instead, let us pick T = at the 
initial slice at i = 0. We then obtain a much simpler set of evolution equations: 



dS^ 

dt 
dT 

'dt 
dKi 

"df 



r 

= K2 



~ e 



2v-2\ 



dr 

2A _2 



1 , 

-ru' + 1 
2 



dS^ 
dr 



(3.9a) 
(3.9b) 

(3.9c) 



+ 47re^^ {l-Cl)p 
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dK2 
dt 



2u-2X 



2 dSi _ dS3 
r dr dr 



+ 2A' - 2u' 



r 



and the Hamiltonian constraint reduces to 



dS3 
dr 



+ 2 



3 • 



STce^" (1 - C,2) p , (3.9d) 



(3.10) 



Note that we still have an equation for T. Our goal is to use the remaining gauge freedom in 
such a way that we can get rid of this equation. This means that we somehow have to make 
the extrinsic curvature variable K2 vanish. To do so, we first solve the Hamiltonian constraint 
( P.IOD for S'^ and plug it in into ( |3.9d| ), which then reads: 



dK2 
dt 



'2v-2\ 



(3.11) 



We finally fix the gauge completely by choosing the lapse 5*1 such that we make the righthand 



side of ( P.llD vanish. This can be accomplished by requiring that S\ should satisfy the following 
ordinary differential equation: 



dr 



^1^' + ^] S3 + Anre'^^C, p . 



(3.12) 



This condition can also be used to replace S[ and SI in ([3.9c[). Moreover, if we use the Hamil 



tonian constraint ( |3.10[ ) to replace the fluid variable p, we end up with just two equations for 
the metric variable S3 and the extrinsic curvature variable Ki. We can even go one step further 
and combine them to yield a single wave equation for S3. However, we do not write down this 
equation because we will present a somewhat simpler system below. 

Before doing so, we would like to prove that radial gravitational waves do not exist. The 
particular feature of this wave equation for ^3 is that its propagation speed is not the speed of 
light alone, but it is multiplied with the square of the sound speed C^. Hence, outside the star 
this equation loses its wavelike character because Cg = 0. In fact, in the exterior region, where 
in addition to = we have A' = —u', the equation for S3 reduces to 



dt^ 



dS3 f , I 

dr \ r 



S, 



(3.13) 



which by virtue of the Hamiltonian constraint ( |3.10| ) with p = becomes 



d'S3 

dt^ 



. 



(3.14) 



This shows us that radial oscillations of neutron stars cannot give rise to any gravitational waves. 
By the way, we could have inferred this much more easily by looking at the momentum 



constraint (|3.7D, which without K2 just reads 



Ki = Aure^" {p + t 



u 



^2v-2\ 



(A' + v')u. 



(3.15) 
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Since u vanishes outside the star, so must Ki and therefore S^. 

To obtain a simpler set of equations than the above mentioned wave equation, we look at the 
so far neglected equations which follow from the conservation law D^T^" = 0. Here, instead 
of p we rather use 



H{t,r) 



p + e 



p{t,r) 



(3.16) 



because the resulting equations are somewhat nicer. We obtain the two matter equations 



'dt ~ ^ ^^d^^^ 
du _ dH dSi 
dt Or dr 



r 



V u 



(3.17a) 
(3.17b) 



Here, we have already made use of T = 5*2 = = 0. We now use the momentum constraint 
( [3.15D to eliminate Ki from equation ( |3.17aD , and our gauge condition ( |3.12| ) together with 
( [3.16D will serve to replace Si in ( |3.17b| ). The resulting equations still contain the metric variable 
53, and in order to obtain a closed system of equations, we need the evolution equation for S^, 
which is given by ( [3.9aD , but with Ki replaced by ( |3.15D . We finally end up with the following 
quite simple set of equations: 



dH 
'dt 
du 
'dt 
dSs 
dt 



2u-2X 



2 9u 



V 



2A' + - 

r 



V u 



87r(p + e)e^''u . 



(3.18a) 
(3.18b) 
(3.18c) 



Additionally, H and 5*3 have to satisfy the Hamiltonian constraint ( P.IOD 

dS3 



dr 



+ 2 



A' ^3 



(3.19) 



In order to obtain physical solutions of the equations, we have to impose boundary conditions at 
the origin and at the stellar surface. At the origin we have to require the perturbation variables 
to be regular. By inspection we find that 



H{t,r) = H%t) + Oir^) 
u{t,r) = u\t)r + 0{r^) 
S3{t,r) = Sl{t)r + 0{r') 



(3.20) 
(3.21) 
(3.22) 



The boundary condition at the surface is given by the requirement that the Lagrangian pressure 
perturbation Ap has to vanish. This condition requires the concept of the radial displacement 
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function ^(t, r), which describes the displacement of a fluid element from its equilibrium posi- 
tion as a function of time and position and will be discussed in more detail in chapter 3. For the 
radial case the relation between ^ and 6ur is given by equation (26.6) of MTW [0: 

% = e^-^^Sur. (3.23) 
ot 

The authors then introduce the renormalized displacement function 

C := r^e-'^e , (3-24) 
for which the Lagrangian pressure perturbation Ap can be written as 

r^Ap = -Tipe-"^ , (3.25) 
or 

where Tip = (p + e) for barotropic equations of state. The condition of vanishing La- 
grangian pressure perturbation at the stellar surface r = R then translates to 

|^(i?) = , (3.26) 

in the case that Tip does not vanish at the surface. For vanishing Tip this is not necessarily true 
anymore, instead, we just have to require the boundedness of ^ and itself. For a more detailed 
discussion of the boundary condition at the surface, see [|TT1]. 

As it is only for the special cases of polytropic equations of state that Tip = at the surface, 
we will always use C'(-^) = as our boundary condition, regardless of the actual equation of 
state. To use this boundary condition for our system of evolution equations ( |3.18| ), we have to 
translate it into a condition for u, which can be easily obtained by making use of ( |3.23| ) and 



(3.27) 



Explicitly we have 



u'{R) = (2X'{R) - u\R) - -^^ u{R) , (3.28) 

which can be used in the numerical code to update the value of u at the surface. This is the only 
relevant boundary condition because the values of the remaining quantities H and 5*3 directly 
follow from the evolution equations. From ( |3.18c[ ) we deduce that S^{R) = 0, and because 
of Cs{R) = equation ( |3.18aD reduces to an ordinary differential equation for H at the stellar 



surface. 

Finally, we should mention that our system is equivalent to equation (26. 19) of MTW, which 
is a single wave equation for the renormalized displacement function (. It can be written in a 



3.2. The eigenvalue problem 
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very compact form with the righthand side having the form of a self adjoint differential operator: 
W^ = —\P^\^QC (3.29a) 



dt'^ dr \ dr _ 
with 

r'^W = (p + e)e^^+^ (3.29b) 

r^P = (p + e) C^, e^^^" (3.29c) 

r^Q = e""^^" (p + e) (^{u'f + 4^ - Svre'v) . (3.29d) 

3.2 The eigenvalue problem 

By applying the Fourier transformation to the numerical evolution we should be able to find the 
frequencies of the eigenmodes. However, it is certainly reasonable to calculate the eigenfre- 
quencies directly from equation ( P.29D , too, by making the harmonic time ansatz 



C(t,r) = e^-*x(r) • (3.30) 

This then gives us a linear ordinary differential equation for x(r) with cu"^, the square of the 
oscillation frequency, as a free parameter: 

Together with boundary condition (|3.26D this defines a Sturm-Liouville eigenvalue problem. 



which has solutions only for a countable set of real eigenvalues cu'^. For tu^ positive, to itself is 
real and thus, the solution is purely oscillatory. However, for negative cu'^ we have an imaginary 
frequency u, which corresponds to an exponentially growing or damped solution. Since the 
general solution is always a superposition of both the growing and the damped modes, this 
means that the occurrence of a negative value of tu^ corresponds to an instability with respect to 
radial oscillations of the stellar model under consideration. For neutron stars this will, indeed, 
happen for central densities eo larger than the critical central density ecru at which the stellar 
mass M as a function of eo has its maximum. In this case the star will ultimately collapse to 
a black hole. For eo = ecrit there must be a neutral mode with the corresponding eigenvalue 

To solve the eigenvalue equation ( p.31| ) numerically, we will write it as a system of two first 
order equations in x and i] := Px''- 

^ = -{uj^W + Q)x- (3.32b) 
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By inspection we find that close to the origin we have x(r) = Xoi"^ + 0{r^) and r](r) = 
r]o + (9(r^). From ( |3.32aD it then follows that the leading order coefficients are related by 3xo = 
r/o/Po, where Pq = {p{0) + e(0)) C2(0) e^(o)+3'^(o). Choosing r/o = 1 we obtain xo = l/(3Po), 
which gives us the initial values for the integration. 

To find the eigenvalues, we will choose some arbitrary to and integrate the equations from 
the origin r = outwards to the stellar surface at r = i?, where we have to check whether 
the boundary condition x'(-R) = is satisfied. If so, the chosen cu corresponds to the desired 
eigenfrequency. 

Numerically, the boundary condition x'{R) = will never be fulfilled exactly. However, if 
we consider x'{R^ ^) a function of to we can find that the zeroes of x'{R^ ^) ^re of first order, 
which means that they are associated with a sign change of x'(-R) Numerically it is quite 
easy to look for a sign change in a given interval, and we can quickly locate the exact value of 
CO up to the desired precision, for example by the method of bisection. Of course, there exist 
other methods of computing the eigenvalues and some of them are compiled in [jTT]]. 

To check our numerical eigenvalue code, we compare our modes with results available in 
the literature. An extended survey of the first two radial modes for a wide range of different 
equations of state was given by Glass & Lindblom [IRp. However, it seems that an error sneaked 
into their computer code, for their equations are correct, but the results are erroneous. Instead, 
we agree with results obtained by Vath & Chanmugam [ [l5ll , who also pointed out that the 
results of Glass & Lindblom are flawed. The strongest argument in favor of their (and therefore 
our) results being correct is that they indeed obtain the neutral modes right at the maximum of 
the mass function A'/(eo). 

In Fig. ^TT] we show the frequencies of the first 5 radial oscillation modes as a function of the 
central density eo. In addition, we include the values obtained by Vath & Chanmugam, which 
perfectly agree with our results. The stellar models were obtained using a realistic equation of 
state described by the model V of Bethe & Johnson [ |75| ] (EOS D in the list compiled by Arnett 
& Bowers ^]). 



3.3 Numerical results for polytropes 

For the first evolution runs we use a polytropic equation of state with F = 2 and k = 100 km^. 



We discretize the system ( |3.18| ) with a two-step Lax-Wendroff scheme (see e.g. [||]), where we 
first perform a half time step to compute intermediate values and then perform a full time step 
to obtain the values at the next time level. 

We show the evolution of a narrow Gaussian profile in 5*3 for three different stellar models. 
Model 1 has a central energy density of cq = 3 • 10^^ g/cm^, which corresponds to a mass of 
M = 1.21 Mq and radius ofi? = 8.86 km. Model 2 with eo = 5.65- 10^^ g/cm^ is right below 
the stability limit and model 3 with eo = 5.67- lO^^g/cm'^ is above it. For models 1 and 2 we 
expect periodic time evolutions with the signal being a superposition of the various eigenmodes. 
Model 3, which is unstable with respect to radial collapse, should show an exponential growing 
mode. 



Our expectations are fully met by the numerical evolution of equations (3.18). In Fig. 3.2 we 
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5.0-10^'* 1.0-10^^ 1.5-10^^ 2.0-10^^ 2.5-10^^ 3.0-10^^ 



Central density in g/cm^ 

Figure 3.1: The first 5 radial pulsation modes as a function of the central energy density for the 
EOS of model V of Bethe & Johnson. The results of Vath & Chanmugam are represented by 
circles. 



show the time evolution for model 1 with its Fourier transformation in Figure The spectrum 
shows that the chosen initial data excite many of the eigenfrequencies of the stellar model; 
at least 15 modes are clearly present. In the spectrum we also include the eigenfrequencies 
computed by directly solving the eigenvalue problem ( p.32D , which agree perfectly. It is only 
for the higher frequency modes that the peaks in the spectrum are systematically located at 
higher frequencies than the actual eigenfrequencies, which is due to insufficient resolution of the 
evolution. By increasing the number of grid points the peaks converge to the right frequencies. 

As the central density increases, the star approaches its stability limit. At the same time the 
frequency of the lowest mode starts to migrate towards zero. The stability limit itself is charac- 
terized by the presence of an eigenmode with zero frequency. As was already stated above, at 
this point the total mass M as function of the central density eg exhibits a local maximum. In 
Figure we show the time evolution of H for model 2. The evolution does not really look 
different from the one for model 1, but in the signal there should be a very low frequent os- 
cillation, which corresponds to the lowest eigenmodes. The Fourier transformation in Fig. ^3] 
confirms the presence of a very low frequency mode, which for this model has slipped down 
to a frequency of = 173 Hz. The second mode resides at the much higher frequency of 
z/2 = 7580 Hz. We should mention that in order to obtain the spectrum in Fig. |33| , where we 
have a resolution of about 10 Hz, we have to evolve up to t = 100 ms. With a time step size that 
is somewhat smaller than 1 ■ lO^^'ms we then need more than one million integration steps! 

Model 3 is unstable, which is nicely confirmed by the evolution. With the eigenvalue code 
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Figure 3.2: Evolution of u for a polytropic stellar model with central density cq = 3-10^^ g/cm 
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.3: The power spectrum of the above wave signal shows that at least 15 eigenmodes are 



Figure 3 

present. The circles represent the eigenfrequencies obtained by solving ( [3. 32] ). 



we find an imaginary frequency with Im(to') = 606 Hz, which corresponds to an e-folding time 
of r = 1.65 ms. In the logarithmic plot of Fig. ^]6| the exponential growth shows up as a linear 
increase in the amplitude. From a fit of an exponential function to the numerical data, we find 
r = 1.69 ms. 



3.4 Getting into trouble: Using realistic equations of state 

So far we have used polytropic equations of state, which are quite decent approximations to 
realistic equations of state as far as general features of neutron stars like mass and radius are 
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Figure 3.4: Section of the evolution of u for a polytropic stellar model with central density 
Co = 5.65-10^^ g/cm^, which is right below the stability limit. 




100 1000 10000 

Frequency in kHz 

Figure 3.5: Fourier transformation of the above wave signal. Here, the lowest mode, which lies 
at vi = 173 Hz is clearly visible. The second mode lies at 1^2 = 7580 Hz. 



concerned. However, it is in particular the oscillations of neutron stars that are very sensitive to 
local changes in the equation of state, which are due to the different behavior of the neutron star 
matter under varying pressure. It is therefore much more interesting to use realistic equations of 
state that take into account the underlying microphysics which determines the state of the matter 
as a function of pressure and temperature. For comprehensive overviews on realistic equations 
of state, see [0,1, |]. 

As was already mentioned in the first chapter, for the sake of simplicity we will resort to 
zero temperature equations of state only. Of course, if we were interested in damping times. 
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which are due to internal friction and other viscous effects that result from the finite temperature 
inside the star, we would have to abandon the zero temperature approximation. 

Realistic equations of state cannot be given in analytic terms over the whole pressure range 
inside the neutron star, hence they usually exist in tabulated form only. To solve the TOV 
equations in this case, one has to interpolate between the given values in order to obtain the 
stellar model with continuous functions of radius r. 

The thus obtained stellar models are not as smooth as those with polytropic equations of 
state, for it is clear that the energy density and other matter functions will have bumps and edges 
as the matter undergoes phase transitions for increasing pressure, which can change its stiffness 
quite abruptly. This is particularly the case close to the stellar surface, where the pressure is 
zero and increases extremely as one moves into the stellar interior. It is clear that it is here 
that the matter undergoes a couple of phase transitions. First the crystal lattice gets destroyed 
and changes into a plasma of nuclei and free electrons. As the pressure further increases, the 
electrons get captured by the protons of the nuclei, which yields more and more neutron rich 
nuclei. Eventually the nuclei start to dissolve and neutrons begin to "drip" out of the nuclei. At 
this neutron drip point, which occurs at 7.8-10^^ dyn/cm^, or at a density of about 4.3-10^^ g/cm^, 
the equation of state has its most drastic change. For increasing pressure the matter then stays 
in the form of a degenerate Fermi gas until it reaches the point beyond nuclear density, where, 
again, phase transitions may take place. 

This is the point where the underlying physics is the least known, and one has to rely on 
physical models with simplifying assumptions in order to theoretically compute the equation of 
state. Of course, different nuclear models lead to different equations of state and it is therefore 
at the nuclear density level and beyond that the equations of state that are given in the literature 
differ the most strongly from each other. For a more detailed description of the physics at high 
densities, we refer the reader to IQ. 
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In the following we will make use of an equation of state called MPA [pO|], which yields 
a maximal mass model of 1.56M0. We can obtain a typical stellar model by taking a central 
density of eo = 4-lO^^g/cm^, which yields a radius and a mass of R = 8.18 km and M = 
1.55Mq, respectively. 

If we try to repeat the evolution of an initial perturbation of the radial velocity field u for 
the above stellar model using the Lax-Wendroff scheme for the discretization of the relevant 
equations ( p.l8| ), we will witness a stellar explosion. That is, after a few oscillations we will 
suddenly find an exponentially growing mode that immediately swamps the whole evolution. 
This cannot be a physical instability, for we have taken a stable background model, and by 
playing around with different resolutions we quickly convince ourselves that it has to be a 
numerical instability, for neither the e-folding time nor the frequency seem to converge to some 
limiting values. 

Apparently, it has to be a problem of the numerical discretization scheme we use, we can 
therefore try to switch to a different one. For instance, we could try to discretize the system 



(3.18) on a staggered mesh. If we write the equations schematically as 



Q = aP' + bP (3.33) 
P = cQ' + dQ , (3.34) 

the discretized form reads 

= Q7 + {PIli/2 - PI1/2) + h~ {PJXy, + Ply,) (3.35) 

p::i/2 = piii/2 + {Qti - Qr') + ^.+1/2^ {q-:^' + Qr') ■ (3.36) 

We can see that Q lives on integer grid points, whereas P lives on half integer ones, that is the 
P-grid is shifted by half a grid point both in r- and t-direction with respect to the Q-grid. In 
( |3.18p we put H on the regular grid and u and S-^ on the shifted grid. 



In Fig. we show the evolution for various resolutions. Again, we find a growing mode 
for low resolutions, but as the number of grid points N is increased, the occurrence of the 
instability gets more and more delayed, and the slope of the exponential growth decreases, too. 
And for = 500 it suddenly disappears. Interestingly, it is even possible to pin the point where 
the instability vanishes down to = 469. For A^ < 469 we still have exponential growth, but 
if we add just another grid point we can evolve the equations arbitrarily long. 

It seems that also for the Lax-Wendroff discretization the instability can be made to vanish 
by increasing the resolution, but the required resolution is extremely high. Too high to allow 
for numerical runs within a reasonable time scale. 

We now try a third possible discretization, where we transform (|3.18|) into a single second 



order wave equation. This can only be performed for either uor S3, but not for H, since in this 
case we cannot totally eliminate the remaining variable S3. Here we choose u, but since the 
resulting wave equation is somewhat lengthy, we define another variable w through 

w{t,r) = re''-^''u{t,r) , (3.37) 
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Evolution using the first order system ( P.18D discretized on a staggered grid for 
ranging from N = 100 up to 500 grid points. The key is the same as in the figure 



for which we have a quite transparent wave equation 
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At the stellar surface, we have the following boundary condition 

{rw)\R) = 0. 



r 



w 



(3.38) 



(3.39) 



We discretize ( [3.38D with central differences, which again will be demonstrated for the schematic 
equation 



Q = aQ" + bQ' + cQ . 
A second order discretization scheme is the following 



(3.40) 



= 2Qr + Q7~' + (Ai) 
h 



2Ar 



(Ar)2 



{Q^^, - 2gr + Qti) 



(3.41) 



In Fig. we show the evolution for different resolutions. The plots are quite similar to Fig. 
Here, too, the instability goes away for resolutions of 500 or more grid points. This similarity 
of the behavior between those two discretizations is understandable, since the discretization of 
(|3.18D on the staggered grid is equivalent to the above discretization of the wave equation. 
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Figure 3.8: Evolution using the wave equation ( |3.38| ) discretized with central differences for 
resolutions ranging from N = 100 up to 500 grid points. 



The question now is, why is there an instability at all? Why is it that for polytropic equations 
of state there is no problem at all, whereas by using a realistic equation of state the numerical 
evolution can blow up? Which part of the equations is responsible for this peculiar behavior? 
This is all part of the investigation in the next section. 



3.4.1 The physical problem 

It becomes clear very soon that the culprit has to be the profile of the sound speed inside the 
star. For if we compute the stellar models with a realistic equation of state, but then replace 
the profile of the sound speed with a profile that results from a polytropic equation of state, the 
instability does not occur any more. 

For polytropic equations of state, the sound speed is a smooth, monotonically decreasing 
function of the radius r that is zero only at the surface of the star, whereas for realistic equa- 
tions of state there are regions in the outer layers of the star, where the sound speed can have 
sharp drops. In particular, close to the neutron drip point, the equation of state becomes very 
soft, which results in a drastic decrease of the sound speed. This sharp drop then can lead to 
the observed numerical instabilities due to an interplay between the boundary condition at the 
surface of the star and the small value of the sound speed in this particular region. 

In Fig. we show the square of the sound speed = ^ for a stellar model using eos 
MPA with central density of cq = 4-10^^ g/cm^. It can be clearly seen that at r 8.06 km there 
is a local minimum of the sound speed, where it drops down to ^ ^ 0.0005. For larger r we 
can see a series of much smaller dips, which is an artefact of the numerical spline interpolation 
between the tabulated points. But the dip at r ^ 8.06km is physical and is present for any 
realistic equation of state. 
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Figure 3.9: Square of the sound speed inside the neutron star model using eos MPA (above) 
and a section near the surface (below). 



3.4.2 A toy model 

The wave equations that describe the non-radial metric oscillations, e.g. the Regge- Wheeler 
equation ( |4.13| ) and the Zerilli equation ( [4-.33D in chapter 4 are of the form 



or 



where 



W 
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_d_ 



or \ or 



dr^ 



c(r) 



d 



dr 



(3.42) 



(3.43) 



(3.44) 
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Here c(r) is the local propagation speed, in this case the speed of light, and V represents a po- 
tential. Besides, the equation can contain additional terms from the coupling to other quantities, 
but they are not relevant for our discussion. By choosing another variable $ = c(r)l>, we can 
further change (|3.42|) to 



52$ 



2/ X^'* T.^ 



(3.45) 



which is the standard form of a wave equation with a variable propagation speed c(r) and 
a potential term V. However, this is not possible for the fluid wave equations. If we let \E' 
represent any of the fluid variables u, w or H, which have to be multiplied by appropriate 
factors e'^ or e^, we cannot find a wave equation of the form ( |3.45| ), with c^(r) now replaced by 
the sound speed C^(r). Instead, we always will have an additional term, which is proportional 
to the first derivative ^' and which we cannot make disappear (we suppress the overall factor 



._,2u-2X 



) 



C^Ar 



'—.-v'—m + Vm . 
ov^ or 



(3.46) 



At first glance this does not seem to be bothersome at all, but it is in the case of the sound speed 
becoming almost zero that the problems start to arise because then the term proportional to \E'' 
starts to dominate. 

For the following discussion we will not use ( |3.46| ) but rather the fluid equation ( |4-.30[ ) 
that governs the non-radial oscillations, which are investigated in chapter 4. From numerical 
experiments we find that the problematic terms in ( 14.30D are the first two ones: 



'2v-2\ 



s Q^2 



+ {Cl {2v' - A') 



) 



(3.47) 



The wave equation ( |3.38D that we use for the radial oscillations is somewhat different, but if we 



were to transform ( p.lSj ) into a wave equation for H instead of u, the first two terms would be 



exactly the same as in ( P.47D . To further simplify ( P.47D , we now set 



g2i/-2A 


= 1 




= 


X' 


= 


1 

V 


= a 



which gives 



Q^2 



(3.48) 



Our goal now is to see if we can get the same kinds of instabilities with this simplified version. 
It turns out that as long as we do not require periodic boundary conditions, the actual choice 



of the boundary conditions does not have a great influence on the stability behavior of (3.48). 
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For the sake of simplicity we therefore choose H{0) = H{R) = 0. We could as well choose 
the first derivative of H to vanish at either boundary. Finally, we set = 1 and discretize our 
toy model ( ]3.48[ ) with central differences: 

I ' (3.49) 

To simulate the dip in the sound speed, we assume a constant sound speed with a small Gaussian 
well with depth h and width w located at Tc. 

c2(r) = Co - /ie"('^)' . (3.50) 

To infer the stability behavior of our numerical scheme, we have to look at the eigenvalues 
of the discretized time evolution operator. Let this operator, which advances data from the time 
level n to the next time level n + 1, be denoted by B. If we furthermore collect all values if" 
on the time level n in the vector if", we can formally write 

if"+^ = Bii" . (3.51) 
Consequently, the inverse of B leads to the previous level 

^n~i ^ B-^ii" . (3.52) 
Our numerical scheme ( |3.49D can be written in vector form as 

ffn+l^jjn^l ^ 2Gii" 

^ (B + B-i)ii" = 2Gii'^, 

where G is the matrix that includes the discretized spatial derivatives. Suppose now that ii" is 
an eigenvector of B with eigenvalue b. Hence the eigenvalue of B"^ is 6"^ and we have 

{b + b-^)W = 2Gii" , (3.53) 

that is, ii" is also an eigenvector of G with eigenvalue 

g = ^{b + b-') . (3.54) 
Solving for b yields 



b = g± ^g'-l . (3.55) 



For stability, we have to have \b\ < 1. If we denote the two solutions of ( P.55D by 



b+ = (3.56a) 



b. = g- ^g^-l , (3.56b) 
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we find that 6_ = 6+^. Hence, if < 1 it is > 1 and vice versa, i.e. we cannot have 
the moduli of both eigenvalues smaller than one at the same time. Therefore it is only for 
1 6+ 1 = 1 6-1 = 1 that we can have stability. This means that b has to lie on the complex unit 
circle: 

6 = 6'"^. (3.57) 
It then follows for g that 

g = - (^e'^ + e-''^) = cos(f) . (3.58) 

Thus, g is real and only assumes values in the interval [—1,1]. 

Hence, if we find any eigenvalue g of G that is either complex or whose modulus is larger 
than one, we always have one |&| > 1 and therefore the numerical scheme is unstable. Now, G 
is a tridiagonal matrix with the band given by the triple 

{Gi,i-i. Gii, Gi^i^i) — 

QcV^ - ^/^Ax (2c^ - 1) , 1-cT, lc'f + lfAx{2c^-l)y 

where we have defined 

A necessary condition for stability turns out to be 

c/ < 1 , (3.60) 
which is the usual CFL criterion. 



Of course, there are a lot of parameters to play around with. It would not make sense to try 
out all the possible combinations. Instead, we choose some parameters to be fixed. For example, 

we choose a = f = 1. Indeed, we have found that the onset of instability is independent of 
/ (always assuming / < 1, of course), which means that smaller time steps with fixed grid 
spacing cannot cure the instability. Too large a value of a combined with too low a resolution 
also can give rise to instabilities. The choice of a = 1 ensures that the scheme is stable if c = 1 
throughout the whole domain. 

We now set cq = 1, and Vc = 0.5, that is, the dip is in the middle of the domain. Interestingly, 
for a given resolution, the instabilities may vanish, if we move Vc towards the left boundary, 
whereas a stable set of parameters can get unstable, if we move Tc towards the right boundary. 
That means if an instability occurs at Tc = 0.5, it will necessarily also occur at Vc > 0.5. (In the 
actual case of a neutron star, the dip sits pretty close to the star's surface). Having now fixed the 
parameters /, a, cq and Tc, we are left with h, w and N. 

Our results can be summarized as follows. First, let us choose w to be so small that the 
considered grids will not be able to resolve the Gaussian. Thus, on our grids we have c = 1 



32 



Chapter 3. Radial oscillations of neutron stars 



1.0014 



1.0012 



1.001 



1.0008 



1.0006 



1.0004 



1.0002 




40 



Cmin = 0.005 -+- 
Cmin = 0.004 -X- 
Cmin = 0.003 * 

Cmin = 0.002 E- 

Cmin = 0.001 



100 



120 



Figure 3 



60 80 
Number of grid points 

10: Largest eigenvalue as a function of the number of grid points for different values 



everywhere except at the middle grid point, where we have Cmin = 1 — h. (We then have to take 
an odd number of grid points, for only in this case there is a grid point located at r = rc = 0.5). 
By explicitly computing the eigenvalues of G we find that they are indeed real, but they are 
not necessarily smaller than one. It is rather the case that for a given grid size N there exists a 
certain critical value Ccru, for which at least one eigenvalue will be greater than one. For Cmin 
only slightly larger than Ccrit all eigenvalues are < 1, whereas for Cmm < Ccrit there is always at 
least one eigenvalue > 1. 

The higher the resolution, i.e. the larger and the smaller Ar, the smaller Ccru- Conversely, 
for a given Cmin, there exists a minimum grid size Ncru, above which all the eigenvalues are < 1 
and below which there will always be at least one eigenvalue > 1. And the smaller Cmin, the 
larger the required grid size Ncru to obtain a stable evolution. However, for Cmin = there 
always seems to be one eigenvalue > 1, regardless of the chosen resolution. Thus, in this case, 
the scheme is unconditionally unstable. 

In Fig. |3.10| we show the eigenvalues of G as a function of the grid size for various values 
of Cmin- It can be clearly seen that for smaller values of Cmm it takes a larger grid size to squeeze 
the largest eigenvalue below one. For Cmm = 0.002 it is for > 61 that all eigenvalues are 
< 1, whereas for Cmin = 0.001 we must have > 121. 

In Fig. p.ll| we show the eigenvector corresponding to the largest eigenvalue for the param- 
eters Cmin = 0.004 and A^ = 21. 

Figures p.l2| and |3.13| show the numerical evolution for Cmin = 0.002 for grid sizes of 
A^ = 59 and A^ = 61, respectively. From Fig. [3.10| we see that for A^ = 59 there exists one 
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Figure 3.1 1: Eigenvector corresponding to the eigenvalue larger than one for grid size = 21 
and Cmin = 0.004. Compare with numerical results in Fig. [3.12 . 

eigenvalue that is > 1, whereas for grid sizes of = 61 or greater, all eigenvalues remain < 1. 
As expected, the numerical evolution shows a growing mode for N = 59, whereas for = 61 
the evolution remains bounded. From the shape of the growing solution we find, indeed, that it 
corresponds to the eigenvector of Fig. p.l 1| . 

The influence of the dip width w is rather counterintuitive. The larger w, the wider and 
smoother the dip. If the instability should be due to the discontinuity in c, it should vanish for 
increasing w. The contrary is the case. For a given resolution and a value of Cmm just above Ccru, 
a widening of the width, i.e. an increase in w will eventually result in a transition to instability. 

We would also like to stress that the instability is not invariably related to the dip, but is 
rather due to the interplay of the dip with the prescribed boundary conditions. For if we were to 
choose periodic boundary conditions, all the previously unstable cases suddenly would become 
stable! This is due to the fact that with periodic boundary conditions our domain is not bounded 
any more, in fact, it becomes infinite. 

This peculiar behavior is apparently present for any numerical discretization, however, for 
the Lax-Wendroff scheme the required resolution in order to suppress the instability can be so 
high that it can prevent one from performing evolutions within a reasonable time frame. It is not 
clear to us why for the Lax-Wendroff scheme the instability is much more persistent than for the 
other two schemes we were investigating. But this shall be of no concern, for we will rewrite 
the equations in such a way that they can be integrated in a stable way for any resolution. 
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Figure 3.12: Numerical evolution for Cmin = 0.002 and = 59. Exponential growth of that 
eigenmode whose related eigenvalue is > 1 is observed. Compare the shape of the mode with 
Fig-FTTl 




Figure 3.13: Numerical evolution for Cmm = 0.002 and N = 61. All eigenvalues are < 1, and 
the numerical solution remains bounded. 
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3.4.3 The solution of the problem: Rewrite the fluid equation 

From the previous section it has become clear that the exploding modes are due to the smallness 
of the sound speed combined with insufficient resolution. Increasing the resolution for a given 
profile of the sound speed or increasing the minimum value of the sound speed for a given 
resolution will eventually result in numerical stability. 

However, for the physical equations the profile of the sound speed is fixed for a given stellar 
model and it is only by further increasing the resolution that we can prevent the occurrence of 
the instability. 

Yet, the nature of the instability is such that we only need the high resolution in the small 
region close to the surface of the star, where the sharp drop in the sound speed occurs. If we use 
a uniform grid, we have to use the same resolution throughout the whole domain even in those 
places where it would not be necessary. And in our case, this would be 99% of our domain! 

The natural way out one might first think of is to locally refine the grid and to provide the 
required resolution only in the region where it is really needed. This could be accomplished by 
fixed mesh refinement since this region is determined by the profile of the sound speed, which 
does not change throughout the evolution. However, we then would have to deal with the tran- 
sition from the coarse grid to the fine grid and vice versa, which might be troublesome. Another 
drawback is that for different stellar models we would need a different grid refinement and it 
would be a matter of trial and error to find the appropriate refinements for a stable evolution. 

Yet, there is a better way out. We can try to find a new radial coordinate x which is related 
to the actual radial coordinate r in such a way that an equidistant grid in x would correspond to 
a grid in r that becomes automatically denser in regions where the sound speed assumes small 
values. A simple relation between the grid spacings Ax and Ar that would have the desired 
properties is 

Ar(r) = C,(r)Ax . (3.61) 



An equidistant discretization with a constant grid spacing Ax would correspond to a coarse grid 
in r for large values of Cg, which becomes finer and finer as the sound speed Cg decreases. 
From ( p.61[ ) we can immediately deduce the form of our new coordinate x as a function of 



r. By replacing in (3.61) the A's by differentials we obtain 



dx 
dr 



1 



Csir) 



(3.62) 



or 



x(r) 



dr' 

c7?] 



(3.63) 



As a consequence, the derivatives transform as 
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and 

^^^C^l,^C4C.y^. (3.65) 

From the last relation we see that the thus defined coordinate transformation will transform the 
wave equation in such a way that the propagation speed with respect to the a;-coordinate will be 
one throughout the whole stellar interior. 

Of course, we have to use ( p.62[ ) with some caution, for if Cs = this transformation 
becomes singular. And this is what happens at the stellar surface. If, for instance, the profile of 
the sound speed is given in the form Cg = Co{l — r / R), where we have Cs{R) = 0, we can 
find an analytic expression for x 

4r) = -^logfl-^) , (3.66) 



Co \ R 

which tells us that at the surface r = R we obtain x{R) = oo. In this case the coordinate 
transformation seems to be quite useless since numerically we cannot deal with a grid that 
extends up to infinity. We would have to truncate it somewhere. But from a numerical point 
of view this is not that bad since going to infinity in the x-coordinate would mean to have 
an infinitely fine resolution in the r-coordinate at the stellar surface. But this is numerically 
impossible as well, so truncating the x-coordinate at some point means to define a maximal 
resolution in r at the surface. 

It should be noted that the above transformation is of the same kind as in the definition of 
the tortoise coordinate r^,, which we will briefly encounter in the next chapter. Here, too, the 
definition 

^ = e'-^ (3.67) 

leads to wave equations with propagation speeds of one throughout the whole domain (cf. the 
definitions of the Regge-Wheeler ( |4.13D and Zerilli equation ( |4.33| )). Thus, we may call x a 



hydrodynamical tortoise coordinate. 

To obtain the background data on the x-grid, we have to solve the TOV equations ( [Z.3p with 
respect to x. Since we also need r as a function of x, we simultaneously solve ( p.62[ ), too. The 
transformed set of equations then reads: 

= cJ + 47rre2\^ (3.68a) 



dx \ 2r 

Tx = ^^(^-^ + 4^^^''^) (3.68b) 

dp dv , , ,^ ^„ ^ 

ax dx 
dv 

— = . (3.68d) 

dx 
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Figure 3. 14: The r-coordinate as a function of the x-coordinate for three different stellar models. 



Of course, we also have to rewrite the wave equation ( |3.38| ) in the new coordinate 



+ a 



r 



2— + 



c. 



CI 



dw 
dx 



-a 



1 

— + - 

C., r 



w 



(3.69) 



Here, the subscript x denotes a derivative with respect to x. The last missing thing is the 
transformation of the boundary condition ( |3.39| ). Unfortunately, we cannot transform ( |3.39| ) in 
a straightforward way since at the surface it is Cg = 0, and therefore the transformation of the 
derivative d/dr = (Cs)^^d/dx is not defined. However, the inverse transformation ( [3.64[ ) can 
make sense if we note that if = then any derivative with respect to x has to vanish. This 
is in particular true for w itself. Thus, at the stellar surface r = R, where the sound speed Cs 
vanishes we can impose the following boundary condition for w: 



dw , 
dx 



\x{R) 



. 



(3.70) 



This corresponds to reflection at a loose end on the x-grid. Actually, this boundary condition 
is more general than the old one ( [3.39D since any finite boundary condition for w' on the r- 
grid would translate into ( |3.7(J| ). Conversely, the above boundary condition ( [3. 70] ) ensures the 
finiteness of w' and w at the surface, which is what actually follows from the vanishing of the 
Lagrangian pressure perturbation in ( |3.25| ). We do not even need to know the actual boundary 
condition for w', the equations will automatically lead to the correct one. 



Figure |3.14| shows the r-coordinate as a function of the x-coordinate for three different 
stellar models. From the curves it is clear that an equidistant grid spacing in x will result in 
an equivalent spacing in r that gets very dense towards the surface of the star, which means 
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15: Spectrum of an evolution using the transformed equation ( [3 .691 ) for the MPA eos 
= 200. The circles represent the frequencies of the eigenmodes obtained from (p.32D. 



that this part gets highly resolved. And this is exactly what we need in order to overcome the 
instability and to obtain a decent accuracy. 

In Fig. p.l6| we compare the evolution of w using the wave equation ( p.38D on the r-grid with 
the transformed wave equation ( [3.69[ ) on the x-grid. We show the propagation of a perturbation 
that starts travelling from the origin towards the stellar surface and then gets reflected. In both 
cases we use a resolution of 50 grid points inside the star. For this low resolution a realistic 
stellar model would cause the evolution on the r-grid to be unstable, hence here we use a 
polytropic equation of state. 

In the lower panel of Fig. p.l6| we evolve w on the x-grid, but we plot it as a function of 
ri = r{xi). We can clearly see that the density of grid points increases as one approaches the 
stellar surface. 
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We have cut off the part close to the surface, since here the value of w drastically increases. 
In the part shown the two evolutions look quite alike and it seems that there is no real advantage 
of one over the other. It is only at the surface of the star that the advantage of the better resolution 
comes to light, which is shown in Fig. p.l7| . Here the amplitude of w is much higher because it 
can be much better resolved. 

Finally, we want to demonstrate the effectiveness of our coordinate transformation by evolv- 
ing an initial perturbation for the MPA equation of state. We choose the central energy density 
to be 3-10^^ g/cm^, which yields a stellar mass of M = 1.49M0. The resolution is N = 200 grid 
points, which would be too low for the other cases to yield a stable evolution. Here we do not 
have any problems, the evolution is stable for any chosen resolution. 

In Fig. p.l5| we only show the spectrum that results from the evolution. We can see many 
sharp peaks, which perfectly agree with the corresponding eigenfrequencies computed by solv- 
ing the eigenvalue equations ( [3.32D . 
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Figure 3.16: Evolution of w on the r-grid (upper panel) and on the x-grid (lower panel). In both 
cases we have excised the part close to the surface. 




Figure 3.17: The excised part of Fig. [3.16| . The region close to the surface gets much better 
resolved on the x-grid (lower panel) than on the r-grid (upper panel). 



Chapter 4 

Non-radial oscillations of neutron stars 



Starting from the general expansions of the metric ( [2.27p and extrinsic curvature ( [2.28p , we 
will choose shift and lapse in such a way that we can obtain the famous Regge-Wheeler gauge, 
which was first introduced by Regge and Wheeler [JTTp in the context of a stability analysis of 
black holes. This gauge was also used by Thome et al. [ [TSl , |19|, ^ |T], ^ in the their 
pioneering papers on neutron star oscillations and still is widely used by many authors. 

The main focus at that time was to investigate the different kinds of oscillation modes that 
a neutron star possesses. Therefore the time dependence of the equations was assumed to be 
given by e^'^*, which reduces the equations to a set of ordinary differential equations. Those 
had to be solved together with the appropriate boundary conditions. The first system that was 
presented was of fifth order, but it was later discovered that it was erroneous, and instead a fourth 
order system could fully describe the non-radial oscillations of a neutron star. Detweiler 
& Lindblom [O] then used this set of equations (which was actually slightly modified due to 
numerical reasons) to compute the /- and some p-modes modes for stellar models with various 
equations of state. 

Chandrasekhar & Ferrari [^Tj], however, chose a different gauge, the so-called diagonal 
gauge, and were able to describe the oscillations in terms of pure metric perturbations simi- 
lar to the black hole case. This enabled them to treat those oscillations as the scattering of 
gravitational waves at the given background metric. 

However, the resulting system of equations was of fifth order and therefore allowed for an 
additional solution that had to be rejected because it was divergent at the origin. The diagonal 
gauge was already used much earlier by Chandrasekhar to describe the perturbations of black 
holes [Q]. Here, too, he obtained systems of equations whose degrees were higher by one than 
the equations that were obtained in the Regge-Wheeler gauge. Of course, they also allowed for 
spurious solutions, which had to be rejected on physical grounds. 

In 1981, Xanthopoulos [ [TS] , [79| ] then showed that the spurious solution (which was later 
called Xanthopoulos solution) could actually be used to reduce the degree of those systems by 
one. However, it was still not clear why in the diagonal gauge one was always led to equations, 
which were higher by one degree than in the Regge-Wheeler gauge. 

Shortly after the papers of Chandrasekhar & Ferrari, Ipser & Price [ P3| ] demonstrated that 
also in the Regge-Wheeler gauge it is possible to construct a fourth order system that does not 
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use any fluid perturbations. Furthermore, they showed [M that the reason of the diagonal gauge 



yielding a fifth order system instead of a fourth order system was that this gauge still allowed 
for an additional nontrivial gauge transformation, of which nobody has been aware before. 

The Regge-Wheeler gauge in turn is complete in the sense that it does not allow for addi- 
tional gauge transformations. The main advantage of this gauge is that it halves the number of 



relevant metric coefficients in ( |2.27D . Since some of the vanishing metric components belong 
to the spatial part of the metric, it is clear that we have to choose shift and lapse in such a way 
that during the evolution those components remain zero. 

If we make the following ansatz for lapse and shift (we omit the indices / and m): 

^1 = -^e^-2"^3 (4.1a) 

52 = 26" K2 (4.1b) 
Vi = (4.1c) 
V2 = e'Ke , (4. Id) 

the evolution equations for the spatial metric coefficients V3, Ti and T3 reduce to 

ps - (4.2) 

|-fi = -2e^ir4 (4.3) 
ot 

|f3 = 0. (4.4) 

In addition the evolution equation for K4 depends only on Ti and V3. Thus, if we initially set 
\/3 = f 1 = f 3 = £4 = 0, we obtain dVs/dt = dfi/dt = d%/dt = dKjdt = 0, which 
ensures the vanishing of all those coefficients for all times. 

This leaves us with only three dynamic metric variables namely the axial variable V4 and 
the two polar variables 6*3 and T2. A glance at the Hamiltonian constraint reveals that in the 
exterior region, the two polar variables 6*3 and T2 are not even independent from each other. 
This means that in the exterior region, axial and polar perturbations each only have one degree 
of freedom. This is in agreement with the general observation that any gravitational wave can 
always be described by a superposition of two independent polarization states. 



4.1 Axial perturbations 

Due to the spherical symmetry of the background the polar and axial equations are totally de- 
coupled from each other. As was already mentioned, the axial perturbations are characterized 
by just one metric variable V4. Furthermore, they do not couple to the fluid of the neutron 
star since energy density and pressure are scalar quantities, and therefore their perturbations 
belong to the polar class. Still, there could be some motion of the fluid, described by the axial 
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coefficient of the 4-velocity, but the dynamical equation for turns out to be 

d 



dt 



uz = 



(4.5) 



which shows that nonzero M3 would describe at most some stationary fluid motion. 

Nevertheless, the axial perturbations are not uninteresting and have been studied in much 



detail [ |42| , pO| , p2| ]. It has been shown that each neutron star model possesses a characteristic 
spectrum of quasi-normal modes, which depends mainly on the compactness of the model. It 
is clear that those modes cannot be excited be means of fluid perturbations of the neutron star, 
but they can very well be induced by impinging gravitational waves or by the gravitational 
potential of a large mass moving on a close orbit. In this thesis, however, we mainly focus on 
polar perturbations, but for sake of completeness we will also write down the equations for axial 
perturbations. 



Yet, we will use a slightly different expansion from the one in ( [2.27D and ( [2.28D , which 
somewhat simplifies the resulting equations. The nonzero axial components of the metric are 
(summation over all / and m is implied): 



(4.6a) 
(4.6b) 



and the ones of the extrinsic curvature read 



kgg kg 

kd,9 kA 



■sin ^OXira sin eWim 
sin 6 Wim sin 9 Xi^ 



(4.7a) 
(4.7b) 



with Xim and Wim defined in Equations ( |A.4D and ( \A.5\ ) of Appendix A. This gives us the 
following set of evolution equations (we again omit the indices / and m): 



'dt 

dt 
dt 



\ dr V r 
/(/ + 1) -2 

dVi 
dr ' 



v. 



and one constraint equation: 



'-Kq = 167re^(p + e)M3 . 



(4.8a) 

(4.8b) 
(4.8c) 



(4.9) 
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The three evolution equations can be combined to yield a single wave equation for the metric 
perturbation V4. If instead of Va we use 



r 

this wave equation reads 



(4.10) 



+ e'"" Utt (p - e) + 



6m + 



where is the tortoise coordinate defined by 



Q 



(4.11) 



dr 



(4.12) 



In the exterior, we have m = M the total mass and e = p = and ( |4.1 1| ) reduces to the famous 
Regge- Wheeler equation 



+ e 



Q. 



(4.13) 



In this case we can give an analytic relation between and r: 



r + 2M log 



2M 



- 1 



(4.14) 



4.2 Polar perturbations 



Here, too, we will use slightly different expansion coefficients from ( [2.27D , which will give us a 
set of evolution equations that is well suited for numerical treatment. We decompose the metric 
as follows (again, summation over / and m is implied): 



a 



Pr 



T 



Im 



+ rS 



Im 



Y, 



Im 










^ '-plm 





r sm 






2 rplra 



\ 



lin 



(4.15) 
(4.16) 

(4.17) 



and for the extrinsic curvature (symmetric components are denoted by an asterisk) 

1 , 



k 



„2A T^lm _d_ 
-^2 dB 



„2A T/-lm _d_ 

-^2 as 



' * r [K^^ - IKY) I Yi^ 

* Or sin^ B [K^^ - 1K\^) 



(4.18) 
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In addition we have the matter variables 




(4.19) 



r 



{Sue, Su^) 




r ^ \ de ' Oct) 



) 



(4.20) 



(4.21) 



It is only through this particular choice of expansion coefficients that we obtain a system of 
equations that can numerically be integrated in a quite straightforward and - this is the main 
point - a stable way. It is a quite well known and particularly bothersome feature of spherical 
coordinates that the resulting equations usually contain divergent terms, which exactly cancel 
at the origin, or quasi- singular terms, which at the origin reduce to indefinite expressions of 
the kind 0/0. For physical reasons the latter have to be finite. If not intuition then at least 
bitter experience tells us that it can be impossible to obtain a numerical stable evolution scheme 
for the "raw" equations without further manipulations. It is clear that a standard numerical 
discretization scheme cannot fully take into account the cancellations that occur on the analytic 
level, which results in severe instabilities at the origin. 

As a first step to overcome those difficulties one therefore has to use special linear combi- 
nations of the variables to avoid the presence of the divergent terms that have to cancel. That is 
the reason why we use the particular combinations of S and T in the expansions of the lapse a 
and the metric component /i„. Furthermore, the scaling with r (e.g. p'^/r instead of just p'™ in 
( }4-.19D ) has been chosen such that the remaining terms, which become indefinite expressions at 
the origin, are only terms that do not contain any derivatives of the variables themselves. With 
this choice all the coefficients but have the same behavior at the origin, namely 



The leading term of is proportional to r . We will discuss the different ways of writing 
the equations and the associated numerical problems in much more detail in the section on the 
numerical implementation. Besides, for the sake of clarity we will again suppress the indices / 
and m throughout the rest of this chapter. 

Our last step before writing down the equations is to replace Ki by the following quantity 



This is necessary in order to get rid of the last remaining singular terms which might be a threat 
to the numerical evolution. 

In this way, we obtain a system of 5 coupled evolution equations, which are of first order in 
time but still second order in space. There are 2 equations for the metric variables S and T and 



Q""(t,r) 



g7(t)r'+i + gf (t)r'+3 + ... . 




(4.22) 
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the 3 more for the extrinsic curvature variables K, K2 and K^: 
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dK 

'dt 
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.2A 



T 



(4.23a) 



S 



(4.23b) 



dT 
'dt 

dt 



dK2 
dt 



K5 



(4.23c) 



+ 



u' X' e 
- + 3- + 2- 



2A 



1 



— e 



2A 



+ 2 {rv' + rA' - 1) 5 



T 



(4.23d) 



2u-2\ 



- Svre^^ (1 - Cll) p 



(4.23e) 



We could easily convert this system into a first order system in time and space by adding another 
two evolution equations for the first derivative S' and T' . However, the form of the first four 
equations ( 14.23a| ) - ( |4.23d| ), which are independent of K2, suggests to rather convert them into 
two coupled wave equations for S and T 



d^S 



•2v-2\ 



^'^-|-(5 ' A') 

Q^2 
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(4.24) 
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d^T 
'W 



•2v-2\ 



dr"^ dr 



2X 



fv' X' e 
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+ 2 irv' + rA' -1)5 
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2A 
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(4.25) 



+ Sttc^'' (1 - d) p 
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which are equivalent to equations (14) and (15) of Allen et al. [Q] (Their variables F and S 



Allen 



are related to ours as follows: F = T and SAiien = e S). As can be seen, the wave equation 
for S ( |4.24D is totally decoupled from the fluid variable p, which only couples to the metric 
perturbation T in ( [4.25[ ). The equation for K2 is only necessary in the interior region, where it 
couples to the hydrodynamical equations, which follow from energy conservation D„T^'^, and 
are given by 



dp 
di 



dui 1 du2 



+ {p + e 



dr r dr 
+ (?>u' - A' + 
dK2 



Ui 



V 



+ e 



2A 



3 



— + (2 + rX') K,--K~ -K, ) + K, 



U2 



(4.26a) 



dui 
~dt 
du2 
'dt 



C. 



2^ 
dr 



(1 + Cl) + {d)') p-\{p + e) (r^f + f + 2rS^ (4.26b) 



cIp - 2 + ^) (^'^ + 



(4.26c) 



Here, we have defined Ui := (p + e) Ui. By introducing the enthalpy perturbation 

H :-- 



ci_ 

p + e 



P 



(4.27) 



the fluid equations assume a more convenient form: 



dH 
'dt 



ct 



dui 
dr 



dK2 
dr 



+ (2 + rX') K2-—K- -K, 



ru'Ko 



+ e 



2u-2X 



(C^ (2z/' - A') - u') 



dui 
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du2 
~dt 



+ e 

dH 
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2u~2X 



X' w 
2- 

r 



2 e 

r 



2A 



U2 



2dS dT ^ ^ 
- - I r^— + — + 2rS 
dr dr 



H-'-{r'S + T) . 



(4.28a) 

(4.28b) 
(4.28c) 



Interestingly, from ( |4.28b| ) and ( |4.28cp it follows that the coefficients ui and U2 are not indepen- 
dent of each other but rather are related via 



Ui 



du2 
dr 



(4.29) 



The above system ( |4.28| ), too, can be cast into a second order wave equation for H, which is 
equivalent to equation (16) of Allen et al. [B4p (the different signs in the terms containing S and 



4.2. Polar perturbations 



49 



T are correct): 



2u~2X 



+ {C', {2u' - A') 



+ (Cr(- + 4- 



-,2A 



+ c: 



7z/ y _ 
2 r r 



,2A 



u') — 

dr 

+ 4 + 7'^ 



^(2r.' + i))(r^5 + T) 



(4.30) 



It thus seems that the polar oscillations of neutron stars can be completely described by three 
wave equations inside the star and two in the exterior region. However, it is even possible to 
further reduce the number of equations, for we have not made use of any of the remaining 
constraint equations. The Hamiltonian constraint relates the fluid variable p to the two metric 
variables S and T: 



+ + 2-2rA' + -e-/(/ + l))5 



,2A 



3^ 

r 



-,2A 



(4.31) 



T . 



It is therefore possible to eliminate p in equation ( |4.25D and thus to obtain a consistent system of 
equations, where both in the interior and exterior the two spacetime variables S and T are used 
to describe the evolution of the oscillations. In the exterior, both S and T propagate with the 
local speed of light e'^"'^, in the interior, however, T then changes its character and propagates 
with the local speed of sound e^'^Cg- 

With p being eliminated, S and T in the interior now become independent variables and the 
Hamiltonian constraint ( |4.31D serves as a definition for p. In the exterior, however, S and T are 
not independent but have to satisfy the Hamiltonian constraint with p set to zero. Unfortunately, 
we cannot use the Hamiltonian constraint to further eliminate one of those variables, but it is 
possible to combine S and T to form a new variable Z (we use the definition (20) of Allen et 

ai. m) 



1 - 



2M 



+ /(/ + !)- 2 + ^ 



2rT' + 



2M-r (2 + /(/ + !)) 
r -2M 



T - 2r^S 



(4.32) 



which then satisfies a single wave equation, the famous Zerilli equation that was first derived in 
1970 by F. Zerilli [El] in the context of black hole oscillations: 



d'^z 



d^z 



Qf2 Qj.2 

Here, we use 2n = l{l 



2e — — --^ Z . 



r^inr + 3M)^ 



(4.33) 



1) — 2, and is again the tortoise coordinate from section 
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The last set of equations that is still missing are the momentum constraints: 

IGne^" {p + e) Ml 



IGvre ^ {p + e) U2 



dr 
+ rK - 

dr 



dr 



- r^K + r(u' + A') Ko - 2K. 
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(4.34a) 



(4.34b) 



They do not provide us with new information since they are equivalent to the time derivative 
of the Hamiltonian constraint in the following sense. Let H, Aii and A^2 denote the righthand 
sides of ( |4-.3ip , ( |4.34a| ) and ( [4-.34bp , respectively. As already pointed out in [P7[], we then find 
that the following relation holds in the exterior region: 



dt^ dr X"^^ r 



2v' 



Ml 



And conversely, we have for the time derivative of the momentum constraints 



^ AW 

dt^^ 



d_ 

di 



Mr 



4r (^"^) 
d 

dr 



M2 . (4.35) 

(4.36a) 
(4.36b) 



This is nothing else but the contracted Bianchi identities for the Ricci tensor and can be checked 
by explicitly differentiating the constraints with respect to t and then making use of the evolution 
equations. 

In the interior, we have already written down the connection between the constraints since 
here the Bianchi identities are equivalent to the conservation of energy-momentum by means of 



the field equations and are exactly given by the fluid equations ( |4.26D or ( [4.28D . 

From the above relations it is clear that if the Hamiltonian constraint is satisfied for all times, 
so are the momentum constraints, and vice versa. 

Let us now turn to the question which system of equations we should use for the numerical 
evolution. The basic idea was to use the first order system that results from the (3+l)-split of 
the field equations. However, because of the instability problems at the origin we had to recast 
the first order system in such a way that it more or less became equivalent to the system of wave 
equations rewritten in first order form. This means that there is no real advantage any more in 
sticking to the first order system, on the contrary, from a computational point of view, it is much 
more efficient to use the wave equations because we need fewer equations. 

If we were to use the first order system Q4-.23| ) together with two auxiliary variables for 
S' and T' and with the three fluid equations ( |4.26D , we would have to solve 10 equations in 
the interior and 6 equations in the exterior. Of course, in the interior we could also use the 
Hamiltonian constraint ( |4.31| ) to eliminate p in ( [4-.25D . In this case, both in the interior and the 
exterior, we would have to evolve 6 equations. 

If we took the wave equations the maximal set would consist of only three equations, namely 
( [4-.24D , ( |4.25D and ( |4.30[ ). Here, too, we can use the Hamiltonian constraint ( |4.31D to eliminate 
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the fluid variable p in ( |4.25| ), which leaves us with two wave equations in both the interior and 
the exterior. Since this is by far the fastest way to evolve the perturbations of neutron stars, we 
will use those equations for the numerical evolution. 

The equation for S, 



2u-2X 



^ + (5^' _ A') — 
g^2 Qj. 
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r J 


/ r 



T 



(4.37a) 



is valid both in the exterior and interior, whereas for T we have to distinguish the two cases. In 
the interior we use ([4-.25[) with p replaced with the Hamiltonian constraint ([4.3 ID 
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(4.37b) 



and in the exterior region we use ( |4.25| ) with p set to zero 



^2u-2X 



+ ( - + 3— 



" \ + {v' - A') 
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dr 



1 



e^^^SL±A\T + 2{ru' + rX' -l)S 



(4.37c) 



If we also used the Hamiltonian constraint ( |4.31D in the exterior, we would obtain an equation 
that has lost its hyperbolic character (just set = in ( |4-.37bD ), which would immediately lead 
to instabilities when numerically integrated. 



In the exterior we could also try to switch to the Zerilli equation ( |4.33D , which would have 
the advantage of being a gauge invariant single wave equation. In addition, for large exterior 



grids this would reduce the computing time by a factor of two. From (4-. 32) we can compute Z 
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from S and T, and it is also possible to invert this expression to yield S and T in terms of Z: 
TMr-2M)Z'.(in<.l)-f (l-^))^ (4.38a, 

N 2 . N \ (4.38b) 



+ 3 — -K^ + l) + --r 3 



where 



6M 

A = /(/ + 1) _ 2 + . (4.38c) 

r 

However, the numerical experiment shows that switching in the exterior of the star from the 
variables S and T to the Zerilli function Z and using the Zerilli equation (|4.33|) to evolve Z 



in the exterior causes a numerical instability. This is because S, T and Z are not just related 
by some linear combination, but the relations involve derivatives and are only valid if S and T 
satisfy the Hamiltonian constraint, which is, of course, not strictly true in the numerical case. 
For instance, if we were to numerically compute Z at grid point i from S and T using formula 
( [4.32D , where we approximate T/ by (Tj+i — Tj_i) / (2Ar) and then in turn compute S and T at 



the same grid point i from Z using formulas ( |4.38D , where again we approximate the derivatives 



of Z with central differences, we would see that the resulting values could differ by quite a 
large amount from the original values we started with. During the evolution, this mismatching 
between those values, which would occur at the point where we switch the equations, would 
rapidly amplify and spoil the whole evolution. 

In pB| ] and [^, Moncrief showed that it is possible to construct two gauge invariant quan- 
tities qi and q2, which completely describe the stellar oscillations inside the star. Moreover, the 
fluid-like quantity q2 vanishes in the exterior region by virtue of the Hamiltonian constraint and 
the quantity Q := qi/A satisfies the Zerilli function. It is even possible to define Q in the stellar 
interior by using the following definition for A: 

A = /(/ + !)- 2e-2^(l + rA') , (4.39) 



which in the exterior agrees with ( |4.38cD . This would mean that in the interior we would have 



one wave equation for the fluid variable qi and another one for Q, which in the exterior would 
automatically transform into the Zerilli equation. This therefore would be the most efficient set 
of equations with the additional advantage of being gauge invariant. Unfortunately, Moncrief 
does not write down the relevant equations, which is quite understandable since in the interior 
they become terribly messy. 

We therefore stick to the above formulation of the perturbation equations as a set of two 
coupled wave equations for both the interior and the exterior region of the star. 
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4.3 Boundary and junction conditions 

There are three boundaries we have to take care of. As was already discussed in the previous 
section, at the origin r = we have to demand all variables to be regular. From Taylor expan- 
sion around r = we then can infer the analytic behavior of the various variables, which is 
proportional to r^"*"^ for both S and T. 

At the outer boundary far away from the star, we require the waves to be purely outgoing. 

The third boundary is the surface of the star at r = R, which is formally defined by the 
vanishing of the total pressure P. Since the perturbations will slightly deform the star, the 
perturbed surface will be displaced by an amount ^* with respect to the unperturbed location at 
r = i?. If the coordinates of the unperturbed surface are denoted by x^j, the vanishing of the 
total pressure P at the displaced surface translates to P(t, x*^ + = 0. Taylor expansion to 
first order then gives 

= p(t,4 + r) 



In the last step, we have made use of the fact that the total pressure P is the sum of the un- 
perturbed pressure p and its Eulerian perturbation 5p. In addition, we have omitted the term 
that contains the product of ^* and Sp since it is of second order in the perturbations. Now, the 
unperturbed pressure p is a function of r only and it is furthermore p{r = R) = 0, hence we 
obtain 



p(x^^) + 5pit,x^j,)+e^pix],). 



(4.40) 



Sp{t, x] 



-Cv{R) ■ 



(4.41) 



Unfortunately, this is not a very convenient boundary condition since we neither use the pressure 
perturbation nor the displacement vector in our set of evolution equations. Therefore we must 
relate this condition to the variables we use. We will try to find a condition that gives us the 
time evolution of Se at the stellar surface. The first step is to use the relation 5e = dp/deSp, 
which gives us 
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6e 



dp dt 
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(4.42) 



The time derivative of can then be related to the r-component of the 4- velocity u^. P3|]: 

-'^^{e''5ur- I3r) . (4.43) 
After expansion in spherical harmonics, we finally obtain 
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(4.44) 
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The equivalent equation for the quantity H as defined in ( [4 .271 ) reads 
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di 
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RK2{t,R) + e 
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(4.45) 



Incidentally, this expression can be derived directly from the evolution equation ( |4.28a| ) just by 
setting to zero. The same is true for the wave equation ( |4.30D . For polytropic equations of 
state it is always = at the surface of the star, hence in ( |4.28a| ) the boundary condition 
( [4.45P is satisfied automatically. For realistic equations of state the sound speed at the surface 
should be that of iron, which is very small compared to the sound speed inside the core, where it 
might reach almost the speed of light for very relativistic stellar models. For practical purposes, 
in those cases we just might as well set (r = R) = 0. 

Let us now turn to the junction conditions at the surface of the star. We will always assume 
that e and C| go to zero when approaching the stellar surface r = i?. If at the surface we had 
a finite energy density (as for example in a constant density model), this would result in dis- 
continuities in A', which in turn would affect the differentiability properties of the perturbation 
quantities. For polytropic equations of state, our assumptions are always fulfilled and for real- 
istic equations of state, both the density and the sound speed are very small compared to their 
values in the core and may thus be confidently set to zero at the surface. 

The continuity of the first and second fundamental forms across the surface ensures the 
continuity of the metric perturbations S, T and S'. The associated extrinsic curvature variables 
K, and K' must be continuous as well. From ( |4.23b| ) it follows that S" is continuous, too. If 
in ( |4.23dD we substitute T" by means of the Hamiltonian constraint ( |4.3ip , we can see that T' is 
continuous. The continuity of T", however, depends on the value of p. If we let the subscripts 
in and ex represent the values for the interior and the exterior, respectively, we have for T" 
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(4.46) 



As we shall see below, for polytropic equations of state p(R) can either be zero, finite, or even 
infinite, depending on the polytropic index T. This has to do with the behavior of the derivative 
of the background energy density e', which appears in the boundary condition ( [4.44] ). 

For polytropic equations of state, it is clear that at the surface, both p and e vanish. The TOV 
equations ( [2.3b[ ) and ( |2.3cD can be combined to yield 



,2A 



p' = — {p + e) (m + 47rr^p) 



(4.47) 



from which we see that p' vanishes at the surface. However, e' can behave in quite different 
ways. From (|2.3c|) we find 
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Close to the surface, it is p ^ e and therefore 
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For a polytropic equation of state the square of the sound speed ^ is given by 

^ = KTe'^-' , (4.50) 
de 



hence 



and 



e' = -u'^ . (4.52) 



Approaching the surface, e ^ and u' will become a constant. Now, from ( |4.52| ) we see that 



the behavior of e' critically depends on the value of the polytropic index T. We can distinguish 
three different cases. For F < 2, we have e' ^ 0, for F = 2 we have e' const., whereas for 
F > 2 we have e' — oo! 

This is somewhat disturbing since for the boundary condition ( |4.44| ) this would mean that 
IpI — > oo, unless the expression in brackets vanishes. Unfortunately, this is not automatically 
guaranteed! Interestingly, the boundary condition ( |4-.45p for H is harmless for all values of F 



since u' is always bounded. But, of course, if we were to compute p from H using 

p = + 'h ^ '^H (4.53) 

we would obtain an infinite value when F > 2 unless H vanishes at the surface. However, as in 
( I4.44D , ( |4.45| ) does not guarantee the vanishing of H, even if H is initially set to zero. Therefore, 
we must ask ourselves what really happens in the case F > 2, where it seems that \p\ oo. 



First, we would like to refer to a paper of Moncrief [p9|], where he discusses a sufficient 
stability condition for the non-radial stellar oscillations. For a polytrope he finds that for the 
potential energy in the vicinity of the surface to be positive it must hold that 

/(/ + 1) - 2 - 47rr2— - > 0. (4.54) 
kF 

Moncrief shows that for / > 3 this condition is always satisfied for 6/5 < F < 2. For Z = 2 he 
obtains 6/5 < F < 4/3. The reason for the lower limit of F = 6/5 is that for smaller F there 
are no bounded stellar models. For F > 2 the above condition will always be violated since the 
last term then goes to negative infinity. 

However, the above condition is only a sufficient condition. But it cannot be a necessary 
condition for stability, for in that case all stellar models with F > 2 would be unstable with 
respect to non-radial oscillations. For / = 2 even the models with F > 4/3 would be unstable, 
which is not the case as mode calculations pSj ] and the direct evolution of the perturbation 
equations show. 
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A possibility to clarify this weird behavior of p is to look at the Lagrangian description of 
the perturbations. By definition, Lagrangian perturbations are changes that are measured by an 
observer who moves with the fluid. Hence, she would compare e.g. the energy density at the 
displaced location + to the original unperturbed value at x'. In mathematical terms the 
Lagrangian energy perturbation reads 

Ae(t, x') = £(t, x' + ii{t, x')) - e{x') . (4.55) 

Here, £ denotes the total energy density. A similar expression holds for the Lagrangian pressure 
change. To compute Ae at the surface, we set = and obtain 

Ae{t,R') = S{t,R' + C) -e{R') 

- S{t,R' + C), (4.56) 

since e(i?*) = for a polytropic equation of state. Furthermore, as we already know that the 
total pressure P{t,R^ + has to vanish and since for the polytropic equation of state it is 
£{P = 0) = 0, we obtain Ae(t, i?') = 0. 

Hence, the Lagrangian energy density perturbation always vanishes at the surface, regardless 
of the actual value of F. 

What then happens to the Eulerian density perturbation? By definition the Eulerian density 
perturbation 5e is the difference between the perturbed energy density £ and the background 
density e at the same location 

5e{t,x') = £{t,x')-e{x') . (4.57) 

It is through Taylor expansion to linear order that we obtain the connection between the La- 
grangian and Eulerian perturbations: 

Ae(i,xO = £{t,x')+C{t,x')-^£{t,x^)-e{x^) 

= 5e{t,x')+C{t,x')^^£{t,x') 

= 5e{t, x') + i% x') A + 5e{t, x')) 

= de{t,x') + C{t,xy{r) + e{t,x')-^^de{t,x') . (4.58) 

The last term usually can be dropped with the argument that it is a product ot the two infinites- 
imal quantities and the gradient of Se and therefore of second order in the perturbations. 
However, in the F > 2-case this argument breaks down at the surface, where the gradient of the 
background energy density e' becomes infinite. If we were to drop the second order term, it is 
clear that the Eulerian perturbation Se would have to become infinite, too, in order to compen- 
sate for the blow up of e' and to yield a vanishing Lagrangian perturbation Ae. 

This shows us that in the F > 2-case the physically meaningful quantity is the Lagrangian 
energy perturbation Ae, which remains bounded everywhere and not the Eulerian energy per- 
turbation Se. However, in our case it is not possible to switch from the Eulerian description to 
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the Lagrangian because the latter is only defined in the stellar interior, where we can define a 
fluid displacement vector. Outside the star we have vacuum, which cannot be displaced, hence 
to describe the metric perturbations, we have to rely on the Eulerian description. Therefore, at 
the stellar surface we would have to switch from the Lagrangian to the Eulerian description. 
But it is right there that the Eulerian concept is misbehaved for F > 2, and we would run into 
the same troubles. 

The whole discussion seems somewhat irrelevant, for in the actual set of equations we do 
not use the Eulerian perturbation of the energy density at all. We have got rid of it by using 
the Hamiltonian constraint ( [4-.31[ ). Also, the gradient of the background energy density does 
not appear anywhere in the equations. However, we have to compute the second derivative of 
T, which, as can be seen from the continuity analysis at the beginning of this section, depends 
on the behavior of the Eulerian energy perturbation p. From ( j4-.46[ ) we infer that in the F > 2- 
case we must have a blow up of T" at the surface. Of course, this is very troublesome for the 
numerical discretization, and even for F = 2 we still have a discontinuity in T", which will 
spoil the second order convergence of the numerical discretization scheme. 

The numerical evolutions indeed confirm the above analysis. By computing p with the aid of 
the Hamiltonian constraint ( |4.3 1| ) we find that for polytropic stellar models with F > 2, p tends 
to blow up at the stellar surface, even if it was there initially set to zero. In Fig. ^]T] we have 
plotted the values of p after a certain time of evolution for three different polytropic indices, 
namely F = 1.8, F = 2.0 and F = 2.2. The corresponding values of k are 0.184 km^ ^ 100 km^ 
and 49600 km^ "^ and were chosen in such a way that for the same central density, we obtain 
models with the same radius. As is evident from Fig. |4TT| , for F = 1.8, p vanishes at the stellar 
surface, for F = 2.0 it assumes a constant value and for F = 2.2 it diverges, which is consistent 
with the foregoing discussion. It should be noted that in those numerical simulations we have 
used a grid size of 6400 points inside the star in order to have a decent resolution. 



4.4 Numerical implementation 

It is not quite straightforward to implement a numerical discretization scheme for the above set 
of equations ([4.37[), and one has to worry about some severe problems that will arise when one is 



not doing the right thing. First, we have to treat the boundaries in a correct way. In the previous 
section we have seen that there are in fact three boundaries one has to deal with, namely the 
origin r = 0, the surface of the star and the outer boundary of the grid. The latter is the easiest 
to handle and could in principle (if computing time does not matter) be totally ignored. 

At the outer boundary we impose outgoing radiation condition, which can be realized in 
the following way: Far away from the star, we know that the asymptotic solution ^(t,r) is 
an outgoing wave with a propagation speed c = e"''^ and an amplitude that scales with some 
power of r: 



\l/(t,r) = r"$(cit - r) 



(4.59) 
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Figure 4.1: Snapshots of the values of p after a certain evolution time for polytropic stellar 
models with three different polytropic indices T = 1.8, 2.0 and 2.2. It is obvious that for 
r = 2.2 p diverges at the surface of the star. 



Here \& stands for either S or T. Working out the different derivatives 
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(4.60) 
(4.61) 

(4.62) 



Knowing this relation can help us to compute ^{t, r) at the boundary grid point at the next 
time level n + 1. If we denote the time step by At and the spatial grid spacing by Ar, we can 



discretize equation ( |4.62D at the intermediate grid point (A^ — 1/2, n + 1/2) with ^ and ^ 
approximated by 
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and solve the resulting equation for the unknown value 
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where 

Ar 

A (4.67, 

If a = this reduces to 

n"-' = n^i + ^ - n) • (4.69) 

This is the case for S, whereas for T we have to set a = 1, i.e. the amplitude of T grows linearly 
with r as the the wave travels outwards. As we shall see this will cause some problems when 
we want to compute the Zerilli function. 

Of course, the approximation of ( |4.62| ) with finite differences will lead to a partial reflection 
of the outgoing wave at the grid boundary. The reflected wave will then travel back inwards and 
contaminate the numerical evolution. The effect of this contamination depends on the choice of 
the initial data and in some cases can be quite interfering. For instance we can choose initial data 
that will cause a quite large outgoing initial pulse of radiation followed by a strong ring-down. 
The amplitude of the final ringing of the neutron star can then be smaller by several orders of 
magnitude. In this case, even if only a very small fraction of the first pulse gets reflected, it 
can strongly affect the later numerical evolution. However, we can also construct initial data 
which produce oscillations with more or less constant amplitude. In this case, the reflected part 
practically does not affect the results at all. 

In any case, to be on the safe side, and if computation time does not matter, one can always 
move the boundary that far afield that any reflections will take too long to travel back to the 
point where the wave signal gets extracted. 

The boundary condition at the surface of the star has been described in much detail in the 
previous section. It has become obvious that for polytropic equations of state with the polytropic 
index F > 2 we have to deal with a discontinuity in T" at the stellar surface. 

This, of course, will affect the convergence of the numerical discretization scheme, because 
for instance a second order scheme converges only in second order if the second derivative is 
continuous. And, indeed, by using a second order discretization scheme that does not take care 
of the discontinuity, we only find first order convergence for F > 2. 

However, from a practical point of view, it does not seem necessary to really worry about 
this fact. First, the use of any polytropic equation of state with a constant polytropic index 
throughout the whole star is unrealistic anyway. If we wanted to obtain more realistic results, 
we would have to resort to realistic tabulated equations of state. And for those, e' and therewith p 
is almost zero at the surface, which results in a continuous T". We therefore discretize equations 
( [4.37D with central differences on the whole domain. 



A somewhat more tricky business is the inner boundary at r = 0. The reason is the fact 
that r = is not a physical boundary but rather a coordinate boundary, which is only due to 
the choice of spherical coordinates, and which is absent, for example, in Cartesian coordinates. 
Hence, at r = we cannot impose physical boundary conditions but we must ask for some 
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regularity conditions the functions have to satisfy. To understand what is going on at the origin, 
we look at a simplified version of equation ( |4.25| ), where we focus only on the troublesome 
parts: 

= ^- — -T . (4.70) 

Obviously the righthand side is regular if and only if T has a Taylor expansion of the kind 

T{t, r) = T^it) r'+i + T[(t) r^+' + ... . (4.71) 
The other valid solution 

T(t,r) = To^(t) r-' + ri^(t) r-'+2 + . . . (4.72) 

is diverging at the origin. Hence, we have to make sure, that during the numerical evolution this 
solution will be suppressed. Now, there are several possibilities to modify ( |4.70| ) by rescaling 
the variable T. For instance, we can introduce T = rT, for which we have the following 
equation: 

= \ — -T . (4.73) 

dt"^ dr^ r dr 

In this case T has to be proportional to at the origin. We also can get rid of the r'-behavior 
by introducing T = r'T. The appropriate equation then reads 

d^f d'^f 2(/ + l)9f 

= h -^^ — . (4.74) 

Here T is finite at the origin and symmetric with respect to the transformation r — — r. Hence 
T'(0) = 0, but (T'/r)(0) is finite, again. However, numerically, we cannot directly compute 
this expression at r = since this would result in 0/0, which, as such, is ill-defined. However, 
we could use the I'Hospital rule to obtain {f'/r){0) = (f "/r')(0) = f "(0). 

If we naively try to discretize any of the above equations, say, with central differences, 
we will positively run into troubles. If we ignore the divergent terms for a moment and only 
discretize the simple wave equation T = T" by means of central differences, we find from 
the von Neumann stability analysis that the Courant number C is one. The Courant number C 
determines the maximal allowed time step size Atmax for a given spatial resolution Ar. Since 
in our case the propagation speed is one, we have the relation C = Atmax/^r. 

As soon as we include the divergent term — /(/ + l)T/r^, this is not true any more and 
instead we find C < 1 to be a monotonically decreasing function of /. For large / this means 
that in order to have stability we have to take very small time steps, which makes the evolutions 
more and more time consuming. 



We can compute the Courant condition in the same way as we did in the section |3^, where 
we discussed the influence of the dip in the sound speed on the stability behavior of the dis- 
cretized fluid equation. We have to compute the eigenvalues of the matrix G, which acts on the 
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4.2: The Courant number C as a function of / for the discretized version of equations 
(P773D (solid line) and (^77^ (dashed line). 



vector T" of grid values at the time level n. Those eigenvalues then have to be smaller than one 
in their moduli in order to insure stability. Here, the only free parameter of G are / and the ratio 
At/Ar. 

In Fig. we show the Courant number C as a function of / for the discretized versions of 
( P30| ) and ^J^) with the following boundary conditions at r = 0: T(0) = 0,f (0) = and 

f (0) = {21 + 3)f "(0). It turns out that C{1) is the same for the versions ( pTTS] ) and ( pJOD . As 
I is increased, the allowed time step size drastically decreases. For / = 4 the Courant number 
C has shrunk by a factor of 2 compared to its value for / = 1. It is interesting to note that for 
small values of / version 04.74] ) allows for a smaller time step size than the other versions, but 
for large I things get reversed. This is probably due to the fact that in ( [4.74D the singular term is 
only proportional to I, whereas in ( [4.70D and ( |4.74D it is proportional to P. 

This reduction of the maximal allowed time step size At^ax is a quite undesirable feature, 
but we could live with it if we only used small values of /. But why wait longer for the numerical 
results than really necessary if there is a quite simple trick that allows us to retain a Courant 
number of, say, C = 0.9 for all values of /? 

This trick consists in moving the boundary condition from ro = to some inner grid point 
Tj = iAr, depending on the value of /. The larger / the higher the number of the boundary grid 
point. Specifically, for a given / we choose the boundary to be at ri^i = (/ — l)Ar. Hence, 
for the numerical computation we ignore all grid points rj with i < I — 1 and set T(r/_i) = 0. 
What this basically amounts to is to cut off the bad influence of the /(/ + l)/r^-term close to the 
origin. Of course, by doing so, we artificially introduce some additional numerical error at the 
boundary, but the numerical experiment shows that with the above prescription this error will 
remain bounded and localized only at the boundary. It thus does not have a bad influence on the 
evolution in the remaining computational domain. 

In Fig. 0| we show the Courant number C as a function of / for different boundary grid 
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Figure 4.3: The Courant number C as a function of / for different boundary grid points. 



points. It can be clearly seen that for fixed I the Courant number C increases as the boundary is 
moved to a higher grid point number. Eventually, C will become larger than our desired value 
of 0.9. 

To see whether this is indeed true for a real numerical evolution, we evolve an analytic 
solution to ( |4.70D with periodic time dependence: T(t, r) = To(r) cost, where To(r) can be 
expressed in terms of spherical Bessel functions. 

In Fig. we show the numerical evolution of T with / = 3 for the different left boundary 
points i = 0,i = 1 and i = 2 and C = 0.9. From Fig. we expect the case where the boundary 
is located at z = 2 to be stable, whereas the other cases should exhibit an instability. This is 
indeed what we find from the numerical evolution. In the upper panel, where the boundary is 
at z = 0, we have a quite drastic growth at the origin. In the middle panel, we have put the 
boundary at i = 1, which reduces the strength of the instability, and finally in the lower panel, 
where i = 2, we have a stable evolution. 

So far we have only discussed the numerical treatment of ( [4-.70[ ). Since the singular term in 
( [4-.73D is the same as in ( |4.70| ), the same procedure is also valid for ( [4-.73D . For practical reasons 
we prefer ( |4.70| ), for we save a first derivative. However, things change a little if we want to 
take equation ( [4-.74[ ). In this case we have a different boundary condition at the origin. Still, it 
is possible to move this boundary condition farther to the right, but even by doing so we always 
have to use a smaller Courant number than in the other cases. Since we would like to have a 
fast numerical evolution, we therefore do not use this version. 

The reader may ask why we are presenting three different versions of the same equation, 
when, at the end, we discard two of them. The reason is that when one derives, for example, the 
equations for neutron star oscillations (or any kind of (linear) evolution equations in spherical 
coordinates), the "raw" form usually is not well suited for the numerical treatment. One then 
has the freedom of rescaling the variables. Our suggested prescription is to choose the variables 
in such a way that they have a r'*^+^^ -behavior at the origin. If possible, the singular terms of the 




Figure 4.4: Evolution of T with the left boundary at different locations. Upper panel: boundary 
at i = 0. Middle panel: at i = 1. Lower panel: i — 2. 
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resulting equations should not contain any derivatives. It then should be possible to cure any 
instability at the origin by moving the boundary from r = so far to the right till the numerical 
scheme eventually becomes stable. 

The whole discussion so far has dealt only with the second order equation. However, it is 
also valid for the equivalent first order system. 



4.5 Convergence and the punishment for violating the con- 
straints 

As was already discussed earlier, any valid initial data have to satisfy the constraint equations. 
Once satisfied at t = the evolution equations then guarantee, by means of the Bianchi identi- 
ties, that the constraints will be conserved for all times. 

Of course, in the numerical evolution of the free evolution system, the constraints will start 
to be violated due to discretization and truncation errors. The degree of violation can then be 
used to monitor the accuracy of the evolution. 

In our case we use the free evolution only outside the star, whereas inside the star we have 
already used up the Hamiltonian constraint to reduce the number of equations. Therefore we 
will monitor the Hamiltonian constraint ( |4.31D only outside the star, where it is p = 0. We do 



not monitor the momentum constraints ( |4.34aD an ( }4.34bD , for we have already seen that they 



are equivalent to the Hamiltonian constraint, and since we use the wave equations, we do not 
explicitly compute the extrinsic curvature variables, anyway. 

The exact evaluation of the Hamiltonian constraint should yield zero, however, due to the 
numerical errors throughout the evolution the righthand side of ( |4.31| ) will deviate from zero. 
If the discretization is consistent and the equations are correct, then the violation of the Hamil- 
tonian constraint should converge to zero as the resolution is increased. With a second order 
discretization scheme the error should decrease by a factor of four if one doubles the resolution. 

To check the convergence of our code, we will monitor the Hamiltonian constraint as a 
function of time at some arbitrary location outside the star. In the following we will use poly- 
tropic stellar models with F = 2 and k = 100 km^. The advantage of this choice is that it 
produces smooth functions for p, e, and Cg, but the disadvantage is the discontinuity in the 
second derivative of T. 

As initial data, let us choose a time symmetric gravitational wave pulse that travels towards 
the star and gets scattered back. In Fig. we first show the evolution of the variable T at 
a fixed location outside the star for the two different resolutions of = 100 and A^ = 800 
grid points inside the star. As the amplitude of the wave signal changes by some orders of 
magnitude, we use a logarithmic scale and show the modulus of T. 

The signal consists of three characteristic features. Because of the time symmetry of the 
initial data the initial pulse splits into one outgoing and one ingoing pulse. Since the obser- 
vation point is situated further outside than the initial pulse, we first see the outgoing part of 
the pulse. This pulse is followed by the scattered and distorted ingoing pulse and finally by 
some oscillatory ringing, which consists of the strongly damped w-modes and the fluid /- and 
p-modes, which have much smaller amplitudes. 
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Figure 4.5: Wave form of T for two runs with = 100 and = 800. The main difference 
shows up in the fluid ringing after t = 1 ms, where there is an obvious phase shift between the 
wave forms that are obtained with the different resolutions. 
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Figure 4.6: Same evolution as in Fig. 4.5 for four different resolutions. In addition to the 



variable T, we show the violation of the Hamiltonian constraint ( |4.3lD . For more details see 
text. 



Up to the point where the fluid ringing starts there is basically no difference in the signal for 
the different resolutions. However, in the fluid ringing there is an obvious phase shift for the 
different resolutions. For low resolutions the frequencies of the fluid modes tend to be larger 
than their actual values. Also, the amplitudes of the higher modes are underestimated in this 
case. 
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Figure 4.7: Evolution for a resolution of = 100 grid points inside the star. The upper panel 
shows S and T, in the lower panel we show the Zerilli function Z. Note that even that S and T 
seem to be quite smooth, the Zerilli function looks quite jagged. 



Let us now turn to the Hamiltonian constraint. In Fig. we show the violation of the 
Hamiltonian constraint for the same evolution as in Fig. P3| , this time for four different reso- 
lutions starting with = 100 and doubling each time. For comparison we also include the 
evolution of T. The evolution of the constraint is similar to the evolution of T, but it also has 
some distinctive features. First of all, the oscillation frequencies are much higher, which is 
what we would expect since the violation of the constraint is mainly due to the residual error in 
the approximation of the derivatives by finite differences, which is proportional to some higher 
derivatives. Secondly, in addition to the outgoing and scattered pulses, there is a third one in 
between, which stems from a partial reflection of the ingoing pulse right at the surface of the 
star. For the very high resolutions the violation of the constraint has a lower limit due to the 
finite machine precision, which manifests itself in the noise that is present for N = 800. 

As one can clearly see, in the logarithmic plot of Fig. the spacing between the curves 
for the different resolutions is equidistant, which means that we can translate each curve to 
match its adjacent one by multiplying it with a certain factor. For the upper part of Fig. W^R this 
factor turns out to be exactly 4, which tells us that, indeed, we have second order convergence. 
However, in the ringing phase, which is shown in the lower panel of Fig. W^, this is not true 



any more. Here we can see that the behavior of the Hamiltonian constraint is quite different 
for the various resolutions. The different curves would not match if we tried to superpose them 
onto each other by translating them in the appropriate way. This is because the oscillation 
frequencies of the fluid modes are slightly different for different resolutions. Of course, they 
start to converge for increasing resolution, but unfortunately the convergence is only of first 
order. As was already explained in the previous section, this fact is due to the stellar surface. 
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Figure 4.8: Same evolution as in Fig. ^77] , this time with a resolution of = 800. Again, we 



show S and T in the upper panel, which appear quite similar to Fig. ^J\. In contrast to the 
jagged appearance of the low resolution Zerilli function in Fig. ^77| the high resolution Zerilli 
function in the lower panel is very smooth. 



10"' 



10"' 



10 



10 ■ 



10' 



10"' 



10' 



10- 



10"' 




Figure 4.9a : Zerilli function evaluated at r 
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where the second derivative of T is discontinuous. A second order discretization scheme can 
only yield second order convergence if the second derivatives are continuous, which is not the 
case. 

Another somewhat nasty feature that goes hand in hand with the violation of the constraint 
is the sensitivity in the computation of the Zerilli function. As was already mentioned in the 
section 4.2, the Zerilli function Z is the single gauge invariant function that fully describes the 
polar perturbations of the Schwarzschild spacetime. It can be constructed from S and T by 
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Figure 4.9b : Zerilli function evaluated at at r = 500 km (upper panel) and r = 1000 km (lower 
panel). 
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means of formula ( |4.32| ), however, it is crucial that S and T satisfy the Hamiltonian constraint. 
Unfortunately, even small violations of the constraint can lead to quite crummy Zerilli functions. 
This can be seen in Figs. and p~8i where we show the evolution of S and T together with the 
resulting Zerilli function Z for the two different resolutions of = 100 and N = 800. Whereas 
the wave forms of S and T are smooth and look almost alike for the different resolutions, the 
Zerilli function Z is quite rough for the low resolution, because of the violation of the constraint. 
For = 800 this violation has decreased enough to produce a smooth Zerilli function. 

Unfortunately, the error in the Zerilli function increases as a function of r. This means that 
even if the Zerilli function computed for N = 800 looks very smooth, it can again become quite 
distorted if we move to larger r. We then would have to further increase the resolution in order 
to obtain accurate enough results. 

This is quite unfortunate since if one is interested in the radiation power of the oscillations, 
one has to compute the Zerilli function at spatial infinity. Of course, this is numerically not 
feasible, hence one would like to extract the Zerilli function as far away from the star as possible. 
But there is the caveat. In the wave zone far away from the star, the amplitude of the variable 
S remains bounded as the wave propagates towards infinity, whereas the amplitude of T grows 
linearly with r. The amplitude of the Zerilli function, of course, has to remain bounded, since 
it is related to the radiated power (see appendix C), which has a finite asymptotic value. It 
is therefore clear that in the computation of Z by means of ( [4.32[ ) the growing behavior of T 
somehow must exactly cancel. But, again, this happens if and only if S and T exactly satisfy 
the Hamiltonian constraint. If there is only a slight violation of the constraint, the cancellation 
cannot occur at a hundred percent level and thus the Zerilli function will eventually start to 
grow for large r. This growth is particularly dominant for the high frequency components in the 
wave form. They get much more amplified than the low frequency modes as we move towards 
infinity. Thus, in a power spectrum of the Zerilli function, which is recorded very far away from 
the star we will have a bias towards the high frequency components, which is a pure numerical 
artefact. 



To illustrate the above discussion, in Figs. 4.9a and 4-. 9b we show the Zerilli function Z 
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Figure 4.10: Zerilli function at r = 1000 km evolved with the Zerilli equation. Compare with 



bottom panel of Fig. [4. 9b 



constructed from S and T at the four different locations of r = 125, 250, 500, and 1000 km for 
a given resolution of 400 grid points inside the star. At r = 125 km, which is the upper panel 
of Fig. |4.9a| , the Zerilli function appears quite smooth. At r = 250 km, the panel below, some 



small bumps in the fluid ringing appear, which get enhanced in the upper panel of Fig. |4.9b 
where Z is extracted at r = 500 km. Finally, in the lower panel of Fig. [4-.9b| panel at r = 
1000 km, we can see that the high frequency components have totally distorted the shape of the 
Zerilli function. 



We have already discussed in section ^]2|that it is not possible to switch from the variables 
S and T to the Zerilli function in the exterior region, because of the numerical instability as- 
sociated with this procedure. The most practical way to obtain a quite reliable Zerilli function, 
therefore, is not to completely switch to the Zerilli function in the exterior region but to addi- 
tionally evolve Z together with S and T. That is, close to the surface of the star we construct Z 
from S and T by means of formula ( |4.32| ) and then use the Zerilli equation ( [4.33[ ) to indepen- 
dently evolve Z parallel to S and T. Of course, this amounts to the additional computational 
expenditure of evolving an extra wave equation, but we get rewarded by obtaining much nicer 
results. 



In Fig. |4.1C| we show the Zerilli function evolved with the Zerilli equation and then extracted 

which 



at r = 1000 km. The difference to the Zerilli function in the lower panel of Fig. |4.9b 



should be the same, is quite striking but plausible. There is no amplification of any high fre- 
quency components, and thus the evolved Zerilli function in Fig. |4.10| remains smooth, even at 
large radii. 

In the following section, where we present the results for various stellar models and initial 
data, we always evolve the Zerilli function together with S and T. 
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4.6 Results for polytropic equations of state 

If not mentioned otherwise, throughout this section we always use polytropic models with F = 2 
andK = 100 km^. 

Before presenting any results we have to discuss the various ways of constructing initial 
data, which boils down to the question of how to solve the constraints. Basically, there are only 
two different classes of initial conditions we can construct. 

First, we could excite the neutron star by hitting it with a gravitational wave. This means 
that in the star we initially set all the perturbation variables to zero. In the exterior we can at 
some finite domain arbitrarily prescribe one of the metric perturbations and use the Hamiltonian 
constraint to solve for the remaining one. For the sake of simplicity, we prescribe T to be 
a narrow Gaussian located sufficiently far outside the star. The Hamiltonian constraint then 
reduces to a first order ordinary differential equation for S. Or, we could prescribe the Zerilli 
function and use (t4-.38|) to compute the appropriate values of 5* and T. In this case we would not 



have to solve any differential equation since by construction, S and T satisfy the Hamiltonian 
constraint. 

The other type of initial data is an initial perturbation of the fluid with the associated per- 
turbation of the metric. This means that we arbitrarily prescribe p and use the Hamiltonian 
constraint to solve for the metric variables S and T. However, this cannot be done uniquely, 
since we have two variables to solve for but only one equation. The two easiest possibilities are 
to set one variable to zero and solve for the remaining one. As we shall see, it can make a huge 
difference in the resulting wave form, whether for the same p, we choose S to be zero and solve 
for T or, vice versa, set T to zero and solve for S. 

Of course, all the above choices of prescribing initial data are pure ad hoc prescriptions and 
have no astrophysical meaning whatsoever. As a first attempt to construct initial data by means 
of a more physical prescription, Allen et al. [BTl] have tried to apply the very successful method 



of constructing initial data for black hole collisions in the close limit to the case of the final stage 
of two merging neutron stars. Still, from an astrophysical point of view, this procedure, which 
works surprisingly well in the black hole case is somewhat questionable in its applicability for 
neutron stars. 

The construction of physically realistic initial data is still a very open field and involves 
a deep understanding of the physics behind core collapses or neutron star quakes (glitches), 
which one assumes to be the most important mechanisms to produce significant oscillations of 
the neutron star. That is why we have to stick to our ad hoc data. Nevertheless, they still can 
give us important insight into what kind of modes can be excited, and what generic wave forms 
look like. 

Moreover, for the sake of simplicity, we stick to time symmetric initial data. In this case, the 
momentum constraints are trivially satisfied since all extrinsic curvature variables are zero. In 



[ pSQ Andersson et al. extend the close limit approach of Allen et al. to the head-on collision 
of boosted neutron stars. Here the initial data are not time symmetric any more and therefore 
they have to explicitly solve the momentum constraints for the extrinsic curvature variables. 
Unfortunately we cannot use their initial data in a straightforward way since their data do not 



conform to the Regge- Wheeler gauge. From section 4.2 we know that in the Regge- Wheeler 
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Polytropic stellar models (F = 


2,K = 100 km^) 


iVlULlCl 


eo [g/cm J 


IVl [IVIqI 


It Jvlll 


R/M 


1 


S.O-lQi'i 


0.495 


11.58 


15.843 


2 


l.O-lQis 


0.802 


10.81 


9.137 


3 


2.0-1015 


1.126 


9.673 


5.819 


4 


3.0-1015 


1.266 


8.862 


4.740 


5 


5.0-1015 


1.348 


7.787 


3.912 


6 


7.0-1015 


1.343 


7.120 


3.585 


7 


l.O-lOi^ 


1.300 


6.466 


3.369 


8 


3.0-101^ 


1.097 


5.211 


3.217 


9 


5.0-101^ 


1.031 


4.992 


3.280 



Table 4.1: List of the used polytropic stellar models and their physical parameters. 



gauge the extrinsic curvature variable has to vanish, otherwise Ti, which has to vanish, too, 
would assume nonzero values. The initial data of Allen et al. have non- vanishing K^. 

Let us now turn to the first choice of initial data, the impinging gravitational wave. In 
previous simulations [Q] it was shown that for a reasonable neutron star model in general a 
gravitational wave will excite the first w-mode, the /-mode and several p-modes. However, the 
authors confined themselves to only one polytropic model with eo = 3-10i5g/cm'^ and consid- 
ered only the quadrupole (/ = 2) case. 

Here we would like to present a much more exhaustive survey of wave forms for a whole 
series of polytropic stellar models, ranging from very low mass up to ultra-compact models. Of 
course, both the low mass and the ultra-relativistic stellar models are quite unrealistic; especially 
the latter ones, for they are unstable with respect to radial oscillations. Nevertheless, it is quite 
interesting to see the change in the wave forms as one moves along the series of different models. 
In addition to the quadrupole case, we also pay some attention to Z = 3 and / = 4. 

Because S and T are gauge dependent quantities and as such not very meaningful, in the 
plots we therefore always show the Zerilli function Z obtained by the prescription of the last 
section. 

As initial data we choose a narrow Gaussian in T centered around r = 50 km. For the runs 
with / = 2, the resolution is 500 grid points inside the star. For / = 3 we have to increase 
the resolution to 2000 grid points inside the star in order to obtain reasonable results. The 
observation point is located at r = 100 km and the outer boundary has been moved far enough 
to the right in order to prevent any contamination from reaching the observation point before 



the evolution stops. The physical properties of the used stellar models are given in table P~T. 

Figure shows the evolution of the modulus of the Zerilli function for the stellar models 
Ml - M3 for / = 2 on a logarithmic scale. The wave forms of the low mass models Ml and 
M2 are quite similar: After the reflected wave pulse has crossed the observation point around 
t = 0.5 ms, the amplitude of the wave signal drops rapidly by several orders of magnitude until 
the final fluid ringing starts, which is dominated by the /-mode. The less relativistic the model 
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the lower the amplitude of the fluid ringing. Also, the frequency of the /-mode decreases as the 
models become less relativistic. Obviously, in all three cases, there is no w-mode presence at all! 
This is somewhat surprising since mode calculations show that there exist w-modes for those 
models. Instead, the wave forms show a tail-like fall-off that is characteristic for the late time 
behavior of black hole oscillations. For black holes, the wave form consists of the exponentially 
damped quasi-normal modes and an additional ring-down, which obeys a power law, and which 
results from the backscattering of the gravitational waves at the curved background spacetime. 
After the quasi-normal modes have damped away, it is this tail that dominates the late-time 
evolution. 

In the case of those less relativistic stellar models Ml and M2, the damping of the w-modes 
is quite strong and we therefore can see a tail-like fall-off before the fluid ringing starts. 

It is only with model M3 that things start to change. Here we can see some timid oscillations 
before the amplitude again starts dropping down to give way for the final fluid ringing. In the 
wave form of model M4, which is shown in Fig. |4.12| , a strongly damped oscillation, which 
corresponds to the first w-mode, is clearly visible. In models M5 and M6, with the latter being 
right above the stability limit, it is even more pronounced. In Fig. |4-.17| we show the power 
spectrum of the wave form of model M6. Also included are the frequencies of the /-mode, the 
first j9-mode and the first w-mode that were obtained by a mode calculation with a program of 
M. Leins [|49|]. The agreement is excellent. 

From comparing the wave forms of models M4 - M7 we can discern that the w-modes 
get more and more long-lived as the neutron stars become more and more relativistic, which 
is confirmed by explicit mode calculations. Finally, for the ultra-relativistic models M8 and 
M9 in Fig. |4.13| the wave forms again change their shapes quite drastically, which hints at the 
appearance of a new effect. In those cases there exists another family of very long-lived modes, 
which are called trapped modes, and were first found by Chandrasekhar & Ferrari for the 
axial perturbations of ultra-relativistic stars. As was discussed in section ^]l] the axial modes 
can be described both within and outside the star by a single wave equation with a potential 
term V, whose shape depends on the stellar model. For less relativistic models V is dominated 
by the centrifugal term /(/ + l)/r^ and is therefore a monotonically decreasing function of r. 
As the central density increases, the stellar models get more compact and the remaining terms 
become more important. Above a certain central density, V starts to develop a local minimum 
inside the star, which then gives rise to the existence of the new family of trapped modes. For 
our polytropic equation of state this happens at a central density of about 1.1 ■ 10^^ g/cm'^. 

For the polar case we have two coupled equations and it is therefore not possible to write 
down an effective potential. But there also exist polar trapped modes for stellar models with 
central densities above a certain value, which is about the same where the axial trapped modes 
start to appear. 

In Fig. |4.18| we show two power spectra of the wave form of model M9, one taken from 
the time t = 1 ms, where the trapped modes dominate, and one taken from a much later time, 
where they have mostly damped away and the fluid modes dominate. Here, too, the agreement 
with the modes obtained by a direct mode calculation is evident. 

In Fig. [4-.19I we compare the wave forms of model M6 for / = 2, 3, and 4. It is obvious that 
the higher / is, the less energy the impinging gravitational wave can transfer to the fluid. The 
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Figure 4. 1 1 : Z = 2 wave forms for the low mass models Ml - M3. 
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Figure 4.13: Z = 2 wave forms for the ultra-relativistic models M7 - M9. 




Figure 4.14: I — 3 wave forms for the low mass models Ml - M3. 
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Figure 4.15: Z = 3 wave forms for the intermediate mass models M4 - M6. 




Figure 4.16: I — 3 wave forms for the ultra-relativistic models M7 - M9. 
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Frequency in kHz 

Figure 4. 17: Power spectrum for model M6. The /-mode, the first p-mode and the first w-mode 
are clearly present. 




5 10 15 20 25 30 



Frequency in kHz 

Figure 4.18: Power spectra of the evolution for Model M9 with / = 2. In the upper panel the 
Fourier transformation is taken for an early starting time and shows the presence of the modest 
damped w-modes. The lower spectrum is for a later time, where most of the w-modes have 
damped away and the very long-lived p-modes prevail. 
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Figure 4.19: Wave forms for / 
transferee! to the fluid. 



2,3,4 for model M6. The higher / the less energy gets 



amplitudes of the fluid modes decrease by about two orders of magnitude for increasing /. 

There is another interesting feature in the wave forms for the ultra-compact stars. By com- 
paring the three panels of Fig. \\-.l3\ we can see that the reflection of the impinging wave packet 
with the following quasi-normal ringing gets more and more delayed as the star becomes more 
and more compact. This feature is also evident in Fig. |4-.16| , which is the / = 3 case. 

For model M7 the reflection occurs at around t = 0.7 ms, for model M8 after t = 0.8 ms, 
and for model M9 it is at t = 0.9 ms. For all models the wave forms up to the time of about 
t = 0.6 ms are quite similar. But for the ultra-relativistic models a gap starts to develop between 
t = 0.6 ms and the point where the reflected wave packet crosses the observer and the quasi- 
normal ringing starts. This is due to the fact that the more compact the model the slower the 
local (coordinate) speed of light (which is given by e'^~^) inside the star and the longer the time 
the wave packet takes to penetrate the star, get reflected and leave the star again. In the time 
during which the wave packet is "trapped" inside the neutron star, there is an exponentially 
damped oscillation with a characteristic frequency that depends almost only on / and M and 
that starts in all cases considered at about t = 0.6 ms. 

The surprise is that the frequency and damping time of this ring-down do not correspond 
to any of the quasi-normal modes of the particular stellar model. They rather correspond to 
the (complex) quasi-normal frequency of a black hole with the same mass as the star! The 
correspondence is not quite perfect, but the frequencies agree within less than one percent, only 
the damping of the black hole is slightly larger than that of the neutron star. Still, this is a rather 
unexpected result and, to our knowledge, has not been reported before. 

This peculiar feature is also present in the axial case, which is apparent in Fig. |4.20| , where 
we evolve the same initial data for two different stellar models. The first stellar model is model 
M9 of table |4?T| , the other one, let us call it M9a, is a polytropic star with F = 2.2 and k = 
28911.95. This choice yields the same total mass for the same central energy density as for 
model M9. The radius, however, is somewhat smaller and given hy R = 4.58 km. In addition. 
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Figure 4.20: Evolution of the same initial data for two neutron star models (M9 and M9a) and 
a black hole with the same mass. 

we also evolve the initial data for a black hole with the same mass as the two stellar models. 

Up to i = 0.6 ms all three wave forms are identical since the exterior is given by the same 
Schwarzschild metric for all three cases. As soon as the wave packet hits the surface of the 
star, things start to change slightly. It is evident that the first ring-down phase is not exactly the 
same in the neutron star and the black hole case. Even for the two different stellar models it is 
not identical. But this should not cause any wonder since we are dealing with quite different 
objects, on the contrary, it is remarkable how close the wave forms are up to the point where 
the wave packet gets reflected in the neutron star case. The most obvious difference is that the 
damping for the black hole is somewhat stronger than for the neutron stars. The frequencies are 
almost the same. 

For a black hole the normalized (2M = 1) complex frequencies of the least damped quasi- 
normal modes for Z = 2 and i = 3 are given both for the polar and axial case as 

1 = 2: wi = 0.74734 + 0.17792i 
1 = 3: wi = 1.19888 + 0.18540i . 

Since the potentials inside the star are somewhat different for the two stellar models and 
quite different from the black hole, it must be the outer parts of the potentials, which are the 

same in all three cases, that are mainly responsible for the first ring down. Concerning the stars 
it is clear that the inner part of the potential is important, too, for it is responsible for the stellar 
quasi-normal modes, especially the trapped modes, which do not exist in the black hole case. 
But it seems that in the black hole case the parameters of the least damped mode, which is the 
most important one, are predominantly determined by the outer parts of the Regge-Wheeler 
potential in the axial case, or the Zerilli potential in the polar case. 

To corroborate our presumption, we look at initial data inside the star, now again for the 
polar case. That is we prescribe a fluid perturbation in p and use the Hamiltonian constraint 
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to solve for the metric variables S and T. As was mentioned earlier, this cannot be done in a 
unique way, hence we consider the two extreme cases. For given p we can either set T = and 
solve for S, or vice versa, set S* = and solve for T. 

Furthermore, we can choose p to have its major contribution close to the stellar surface, that 
is at the peak of the potential, or closer to the stellar core, that is right inside the dip of the 
potential. In Figs. |4.21| and |4.22| we show the initial data and the evolution for both cases. 

It is interesting to see that in the first case the wave forms for the two kinds of initial data, 
where we either solve for 5* or for T, are almost identical. Furthermore, it is evident that there is 
no sign at all of a black-hole-like ring-down phase. The spectral analysis of the resulting wave 
forms shows that there are predominately fluid modes and only some of the first long-lived 
trapped modes present. 

In the second case, things have changed dramatically. Here the wave forms for the two 
kinds of initial data are totally different. Those where we initially set T = and solve for S 
lead to a wave form that is somewhat similar to the case of an impinging gravitational wave. 
The frequency of the ring-down after the first pulse is, indeed, close to the first quasi-normal 
mode of an equal mass black hole. It is only after the second pulse that we find the characteristic 
quasi-normal modes of the neutron star. The spectral analysis reveals that both the trapped and 
fluid modes are excited. 

However, if we now initially set 5 = and solve for T, the resulting wave form is quite 
different. There is no trace of any black hole frequency whatsoever, and in the spectrum we 
scarcely find any of the higher trapped modes. All in all the wave form is quite similar to the 
case where the fluid perturbation is located near the center of the star. 

It should be mentioned that a long term evolution shows that both wave forms start to look 
alike. This can be explained in such a way that for both kinds of initial data the fluid modes 
are excited to the same extent since we prescribe the same initial fluid perturbation p, but in the 
case where S differs from zero and T = 0, we have an additional excitation of the iw-modes, 
which because of their strong damping will quickly fade away. Eventually, both wave forms 
contain only the long-lived modes and therefore will have the same time dependence. 

Thus, it seems that in the interior it is the variable S which is mainly responsible for the 
metric modes, and T is responsible for the fluid modes, even if T is actually a metric perturba- 
tion. However, it has already become clear that in the interior T more or less assumes the role 
of the fluid variable. 

We demonstrate this effect more clearly for a less relativistic stellar model, which does not 
have any trapped modes. Hence, here any excited w-modes will immediately damp away and 
give way to the long-lived fluid modes. In Fig. |4.23| we show the evolution for a given initial 
fluid perturbation p, where we either solve the Hamiltonian constraint for T or for S. In the 
former case the resulting wave form is composed of fluid modes exclusively, whereas in the 
latter case we see a short burst of radiation, which then fades away very quickly. The remaining 
fluid oscillations are then almost identical in both cases. 

It thus seems that the prescription of p determines more or less exclusively how much energy 
we initially put into the fluid modes. In addition, the value of S determines how much energy is 
put into the ui-modes. The remaining variable T is then fixed by the Hamiltonian constraint. If 
we set S* = 0, we do not obtain any u'-modes at all. If we increase S, more and more energy is 
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Figure 4.21: Evolution of the initial data which are located inside the potential. The upper 
panel shows the initial fluid perturbation. In the middle we show the two possible solutions of 
the Hamiltonian constraints for the metric variable S with T = and for T with 5 = 0. The 
lower panel shows the evolution of both kinds of initial data. 




Figure 4.22: Evolution of the initial data which are located outside the dip of the potential. As 
in Fig j4.21| the upper panel shows the initial fluid perturbation. In the middle we show the two 
possible solutions of the Hamiltonian constraints for the metric variable 5* with T = and for 
T with 5 = 0. The lower panel shows the evolution of both kinds of initial data. 
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released through w-modes without significantly changing the energy of the fluid modes because 
of the weak coupling between gravity and matter. 

Of course, we cannot choose arbitrarily high values for 5* without affecting the energy of the 
fluid modes. If we increase S by too large an amount, the coupling will transfer the energy of 
the metric into the fluid. Hence, in this case, the excitation of the fluid modes after the ly-modes 
have radiated away will be much stronger than for small values of S. 

It thus seems that the metric variable S has a quite important influence on the form of the 
gravitational waves that are emitted by neutron stars. If we assume the oscillations of a neutron 
star to be excited by either some star quakes or by the precursing core collapse in the stage of 
the birth of the neutron star, and not by an impinging gravitational wave, then the emitted wave 
form will strongly depend on whether or not this particular physical process was able to produce 
some non-vanishing S. From the decomposition of the spatial metric perturbations ( |4.17D we 
see that T is similar to a conformal factor, for if we set 5 = the perturbed metric ( ]4.17D is 
conformally flat. The value of S is therefore a measure of the deviation from the conformal 
flatness. Hence, if the physical process somehow conserves the conformal flatness of the metric 
to some degree, it is clear that one cannot expect to have a significant excitation of w-modes. 

Andersson & Kokkotas [ |55| ] have shown that by extracting the frequencies and damping 
times of the first w-mode and the /-mode in a wave signal one can reveal important stellar 
parameters such as mass and radius, which then can be used to restrict the set of possible 
equations of state. But this procedure stands and falls with the presence of w-modes in the 
signal. It is therefore important to investigate the processes that lead to oscillations of neutron 
stars in greater detail on their ability to excite ly-modes in a significant manner. 

It is also clear that any construction of numerical initial data that rely on conformal flatness 
as in [571, will suppress the presence of w-modes. 
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4.7 Using realistic equations of state 

When we try to switch to a realistic equation of state, we will run into the same numerical 
troubles as in the radial case in chapter 3, i.e. we will have instabilities that are associated with 
the dip in the sound speed at the neutron drip point in conjunction with too low a resolution. 

This is because the structure of the fluid equation ( |4.3U| ) is of the same kind as in the radial 
case. Here we have found a convenient way to get rid of those problems by introducing a new 
"hydrodynamical tortoise coordinate" which stretches those parts of the neutron star where the 
sound speed assumes low values. 

Of course, this transformation is only defined in the stellar interior for the fluid equations, 
for it is only there that the propagation speed is the speed of sound (apart from the factor e'^^'^). 
If we still wanted to use the system of equations ( |4.37D , where in the interior T plays the role 
of the fluid, we would have to switch from the x-grid in the interior to the r-grid in the exterior, 
which is somewhat inconvenient. It is therefore much more natural to explicitly include the 
fluid equation ( [4-.30[ ), which is defined in the interior only. Thence, it is only ( |4.30| ) which will 



be transformed according to (3.62), whereas we keep the wave equations for S and T as they 



are given in ( |4-.24| ) and ( |4.25D . 

However, this means that in the interior we have to simultaneously evolve the fluid variable 
if on a different grid than the metric variables S and T. Because of the coupling at each time 
step we have to interpolate H from the s-grid onto the r-grid in order to update equation ( [4-.25D , 
and, vice versa, both S and T from the r-grid onto the x-grid in order to update equation ( |4-.30| ), 
which can easily be done using spline interpolation. 

The transformed fluid part of the fluid equation ( |4.30D reads 



2u~2X 



dH 

dx 

+ 



ra 



+ 2 



H 



(4.75) 



The subscript x denotes a derivative with respect to x. As in the radial case, the boundary 
condition for H at the surface of the star simply translates into 



dH 







(4.76) 



Apparently, it is the region close to the surface, which is mainly responsible for the fluid modes. 
To obtain the right mode frequencies in the power spectrum, it is important to have high resolu- 
tion close to the surface of the star and this is what can be accomplished by our new coordinate 
X. In Fig. [4-.24| we show the power spectra of wave forms obtained from runs with different 
resolutions. The initial data always were some fluid excitations at the center of the star together 
with 5 = 0. In this case we suppress any w-mode contribution. 

It is apparent that even for the quite moderate resolution of 50 grid points on the x-grid 
inside the star we obtain very accurate frequencies for the first couple of fluid modes. On the 
r-grid, however, a resolution of 200 grid points is still not enough to obtain the same accuracy, 
and the peaks of the higher p-modes in the spectrum are still quite far off their true values. 
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Of course, those results were obtained with a poly tropic equation of state, since otherwise we 
would not have been able to do the evolution on the r-grid because of the occurring instability. 

After all those pleasant features of our new a;-coordinate, we should also mention one draw- 
back that comes along with the transformation of the fluid equation, and which affects the use of 
poly tropic equations of state with F > 2. In this case, we already know from the discussion of 
the boundary condition that the perturbed Eulerian energy density at the stellar surface can ei- 
ther be finite (F = 2) or even infinite (F > 2). This then goes hand in hand with a discontinuity 
or an infinite jump in T". 

In the discretization of the wave equations ( [4-.37P we more or less ignored this fact and 
computed T" everywhere, the stellar surface included, by means of central differences with 
second order accuracy. Of course, for F > 2 the neglect of the appropriate treatment of the 
discontinuity was punished by the reduction of the order of convergence from second down to 
first. But otherwise we did not have any serious disadvantages. 

Unfortunately, with the new coordinate things get much worse. Because of the better reso- 
lution of the fluid at the surface, the discontinuity in T" will be much more pronounced. In our 
case this then leads to some artificial reflections at the surface, which in the case of F > 2 are 
truly severe. In Figs. [4.25| and |4.26| we show the wave forms of S right at the stellar surface for 
different poly tropes with F =1.8, 2.0 and 2.2. Figure \^.25\ shows the evolutions on the r-grid 
and Fig. |4.26| on the x-grid. 

Whereas for the conventional discretization in r the main inaccuracy of low resolutions 
is the phasing, we can see that using the x-coordinate we get some ugly additional bumps in 
the wave function, which start to vanish as the resolution is increased. For F = 1.8, T" is 
continuous at the surface, thus there is only a small reflection for the lowest resolution of 100 
grid points inside the star. For resolutions of 400 or more, the wave form seems to have more or 
less converged to its proper shape. For F = 2, the bumps for low resolutions are much higher 
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Figure 4.26: Same evolutions as in Fig. |4.25| , but this time on the x-grid for T = 1.8, 2.0 and 
2.2, each with different resolutions. Again, for T = 1.8 the convergence is the best, whereas for 
r = 2.2 it is very bad. 
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and it is only for the quite large number of 1600 grid points inside the star that the wave form 
coincides with the one from the r-grid. But for T = 2.2 things seem to be quite hopeless. 
Even for 1600 grid points the shape of the wave form still has almost nothing to do with the 
corresponding one of the r-grid. One would have to double the resolution some more times 
before one could obtain an acceptable wave form. 

This seems to be quite a nasty feature of the new set of equations and it might be worth trying 
to implement an appropriate handling of the discontinuity. However, for our purposes this is not 
necessary. The main motivation for transforming the equations was to get rid of the instability 
that occurs when we switch to realistic equations of state. And this can be accomplished with the 
new version of the equations. Moreover, for realistic equations of state, the perturbation of the 
energy density p is almost zero at the surface, and thus, there is only a negligible discontinuity in 
T". Thus, we will not encounter any of those interfering reflections, and it is possible to obtain 
reasonable wave forms without having to rely on unacceptably high resolutions as it would be 
the case with the old discretization in order to overcome the instability. 

In Fig. [4.27| we demonstrate the usefulness of our coordinate transformation. We show the 
evolution of a sharp peak in the fluid perturbation p, which leads to the excitation of a multitude 
of fluid modes. This is convincingly confirmed by the power spectrum which shows 37 modes 
in the interval up to 100 kHz. And they agree with the modes which are obtained by the direct 
mode computation. 
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Figure 4.27: Upper panel: Evolution of a sharp initial fluid pulse using the MPA equation of 
state. The power spectrum in the middle panel reveals that dozens of p-modes are excited. In 
the lower panel we include the modes computed with explicit mode calculation. 



Chapter 5 

Shooting particles at the neutron star 



In the previous chapter we demonstrated that it is possible to excite the oscillations of a neutron 
star and that the resulting gravitational wave signal contains the expected quasi-normal modes. 
However, all the chosen initial data were purely arbitrary and void of any astrophysical meaning. 
This arbitrariness will now be removed by simulating a physical process that could excite the 
modes of a neutron star. We will use the gravitational field of a moving mass ji to perturb the 
neutron star, hoping that this will lead to excitations of some of its eigenmodes. Since this 
investigation takes place in the framework of perturbation theory, we have to require that the 
mass /i is much smaller than the mass M of the neutron star, or yu/M -C 1. Furthermore, we 
will treat the orbiting mass /i as a point mass. In this particle limit, the gravitational field of the 
mass jj, can be considered as a perturbation of the spherical background field that is due to the 
neutron star. The particle will move on a geodesic in the background metric since deviations 
thereof due to gravitational radiation are effects of second order and will therefore be neglected. 

We will not consider collisions of particles with the neutron star because it is not clear at 
all how to treat the impact and subsequent merge of the particle with the neutron star. Therefore 
we will only focus on circular and scattering orbits. Previous studies of excitations of neutron 
stars by orbiting particles were performed in the frequency domain [ |55| , ^ Since we 
already have the evolution code, we will include the particle and do the explicit time integration. 
As we shall show, considering the dynamics in the time domain has as a consequence that we are 
forced to "smooth out" the particle; that is, the (5-functions in the sources of the equations due to 
the presence of the particle are approximated by Gaussians. We show, however, that our model 
of the particle is self-consistent and convergent. Another important issue is the prescription of 
appropriate initial data that satisfy the constraints. We will face the same problems as in the 
previous chapter, since the constraints do not provide us with a unique way to find valid initial 
data. However, by resorting to the flat space case we can find analytic initial data which can 
serve as a good approximation in the case that the particle is initially far away from the neutron 
star. 

We now proceed how to mathematically describe the particle and how to incorporate it into 
the perturbation equations. 
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5.1 Adding a particle 

The energy-momentum tensor for a point particle with mass fi moving on a geodesic X'^(r) in 
a Schwarzschild metric is given by 



/^T7r^^(^-^W)^(cos^ -cose(t))(5(0-<l>(t)) , (5.1) 



where is the particle's 4-velocity 

f/'^ = , (5.2) 

and T the particle's proper time along its trajectory. We now orient the coordinate system 
in such a way that the particle's orbit coincides with the equatorial plain of the neutron star 
(6 = 7r/2). Besides, the particle's coordinate time T is identical with the time coordinate t of 
the spacetime in which it moves. Therefore we will use t to parametrize the path of the particle 

X\t) = [t, Rit), 7r/2, 

From the geodesic equations = U^DylJ^ = we find: 

^ = e^^E (5.3a) 

CLT 

f^V-^^-e-fl + ^l (5.3b) 

(5.3c) 



_ L 
1^ ~ 



where E and L are the energy and angular momentum per unit mass of the particle, respectively. 
We also recall that for the Schwarzschild metric we have 

9,. 9X 2M 



We can use ( |5.3a| ) to eliminate the proper time r from equations ( |5.3bD and ( |5.3cD : 



(5.5b) 



\dt J V \ 



dt R^E ■ 

On the one hand those equations can be used to replace the quantities vr = ^ and t>$ = ^ 
in the source terms of the particle. But on the other hand we also have to explicitly solve them 
for the particle's trajectory coordinates R(t) and since we need those coordinates in the 
(5-functions 5{r - R{t)) and 5(0 - 
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The perturbation equations for the isolated neutron star are given by ( p.2Up and ( [2.2ip . In 
the vacuum region, the matter terms in ( |2.21| ) vanish, of course. However, when we add the 
particle, we have to include its energy-momentum tensor in ( |2.21D , which then reads 



dthj = -didja + T^dka + ST^dke" 

1 ,.„^ \ (5.6) 
2 



To obtain the initial data, we also have to modify the perturbed constraint equations ( |2.24D and 
( |2.25| ), which then read 

g'^6Rij -h'^Rij = levre-^'^roo (5.7) 
gjk ^Q.j,^^ _ g^j,.^ _ Y\^k,i + r^.A;,,) = Sjre-'Toi . (5.8) 

In writing these constraint equations, we have assumed that initially the fluid of the neutron star 
is unperturbed. 

To simplify the above set of equations, we again want to get rid of the angular dependence by 
expanding the equations in Regge- Wheeler harmonics {[yi^]fj.uid, <P)}a=i,...,w- We then obtain 
equations for the expansion coefficients, which only depend on r and t. For the metric part of 
the equations the decomposition can be done in exactly the same way as in chapter 2, and the 
energy-momentum tensor of the particle can be written as 

oo I 10 

T,^it,r,e,<P) = E E Y.^Ait,r) [ytUiOA) . (5.9) 

/=0 m,=-l A=l 

Unfortunately, the Regge-Wheeler harmonics y'^ do not form an orthonormal set, and it is thus 
not possible to directly obtain the coefficients r) by means of the orthogonality relation. 
We rather have to take some detour. We have to construct an orthonormal set of tensor harmon- 
ics y^, which can be expressed as some linear combination of the Regge-Wheeler harmonics 

y Im' 

10 

yL = Y.^AByEn- (5.10) 

3=1 

We then compute the coefficients r) corresponding to the orthonormal set y^ by using 
the orthogonality relation 



i'T = / (3^zt)*^%.dl]. (5.11) 
Finally, the desired coefficients r) can be obtained by 



10 

tA = Y.^BAit . (5.12) 

B=l 
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which then will be plugged in our equations. 

One convenient orthonormal set can be found in Zerilli [|6T1]. In Appendix A we will list 



the complete set and explicitly show the relation to the Regge-Wheeler harmonics. Li Appendix 
B, we then demonstrate how to derive the coefficients t^X^{t, r). 



If we want to include the particle terms in the equations for the extrinsic curvature (|5.6[), we 
have to use some caution. In chapter 2 we saw that the expansion of ( pT^ into tensor harmonics 
lead to evolution equations for the six coefficients fC-*". By choosing the Regge-Wheeler gauge 
together with the appropriate initial data we could reduce the evolution equation for to the 
trivial case §iKj[^ = 0. This then ensured the vanishing of K^^' for all times. 

However, in the presence of the particle this is not true anymore. Instead of a zero on the 
righthand side, we have a source term 

= -STTtf . (5.13) 

But this means that during the evolution i^l™ will start to differ from zero, which, because 
of (|4-.3[), in turn would make T/™ non-vanishing. This is somewhat unfortunate but it can be 



remedied by choosing a different lapse function a. If we pick 

1 /rplm \ 

a = --eM -- + r5'™ + 167r4'"j , (5.14) 



the last term exactly cancels the righthand side of ( |5.13D and the vanishing of both Ti and 
can be guaranteed. 

Again, we have to distinguish between the axial and the polar perturbations. It is only the 
polar part which will be studied in this chapter, nevertheless, for the sake of completeness and 
because the inclusion of the particle terms in ( p~8| ) - ( |4~9| ) is straightforward, we also present 
the axial equations for the particle: 



dt \ dr \ r 

at 

dKe _ dVi 
dt dr 



Sntio . (5.15c) 



The momentum constraint relates the extrinsic curvature coefficients to the particle's source 
term via 

dKs 2 /(/ + 1) - 2 

— ^ + -^3 - ^ ^ = IQire^H, . (5.16) 
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The relevant source terms of the particle read: 



tiQ = -e 



2u 



d 



(5.17a) 
(5.17b) 
(5.17c) 



For radial infall, it is L = and all source terms vanish. Hence, in this case the radiation is of 
even parity only. 

We can now use the system of equations of chapter 4, and in ( [4.23b| ) and ( [4-.23d[ ) we just 
have to add the following source terms: 



dK 
'dt 



dt 



■ + 1677- 



t'' - 2t' + 5u' 



1 



v' 3 - e- 



2A 



,2A 



^8 ^^9 



■ + IGvre 



J2\ 



- 2^6 + 2 ( z/' - M ^8 + — - ^e^^t 



'-8 i i-g o ^ 

r 2 



(5.18a) 



(5.18b) 



The appropriate source terms are derived in Appendix B and are given by 



^5 
^6 



.6A 



t = -e 



r^El{l + l){l-l){l + 2) 

r^E{l-l){l + 2) 
r'^E 



6{r-R{t))Y, 



Im 



With the explicit form of the source terms, ( |5.18a| ) and ( |5.18bD read: 



dK 
'dt 



2 2y^J' 

iire 



2e^^EvH + im , ( 2r5' + (e^^ - 7)5 

r(n + 1) V 



+ 



r'^Eni 



im'{n + 27) - {n + l){2n + 27)) 5 

(ra + 1) \ ' 



+ (ra + 1 — m 



Y, 



(5.19a) 
(5.19b) 

(5.19c) 

(5.19d) 
(5.19e) 



(5.20) 



Im 
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dt 



r 



E rin + 1 



+ 



e'^'in + 1 - m") {r5' - 65) 



r'^En{n + 1 

+ {{n + l)(2r2 + 3) - m\n + 3)) 5 



Y, 



Im ' 



where we have defined 



n ^ -/(/ + !)- 1. 



For the Zerilli equation ( |4.33D we obtain the following source term: 

= ••■-levre^^— -J^ -le'-'iL' + ry 

r^EA{n + l)[ ^ ' 

f r2 

+ irn^iln + 3) - 4n(n + 1) - 3 + €^''(3 - 2n - 3m^)) 

\2rn ^ 



Imie^^vLE + 12 



A 



r(n+ 1) + 3MU 



1^: 



Zm 



The Hamiltonian constraint ( p77| ) with the particle term reads 



T" + vT' - rS' 



V 



+ e 



2A 



T - ( 2rz/' + 2 + + 1) ) 5 



and the momentum constraints lead to the following two equations 

r/sTg - -e2^/(/ + 1)7^2 - r'^K - (rz/' + 1) 

Here, the appropriate particle terms are 



[lE 



b{r - R{t))Y;^ 



e'^^^^S{r-R{t))Y,l 



(5.21) 



(5.22) 



(5.23) 



(5.24) 



(5.25a) 
(5.25b) 



(5.26a) 
(5.26b) 
(5.26c) 



In the momentum constraints (|5.25[) the quantity K2 still appears, which, again, is not used in 



the evolution equations. Therefore, we will get rid of it by differentiating (5.25a) with respect to 
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r and using ( |5.25a| ) and ( |5.25b| ) to eliminate K2 and K2. The resulting equation is then second 
order in K5 and reads 



K': + u'K' 



V 



„2X 



i{i + r 



K5 -rK' -[ 2rv' + 2 + ]^e^\l + 1) ) ^ 



-Svr + 2 (rz/' + 1) ts - e'^ ^^^^^^ tg^ 



(5.27) 



It should be noted that this equation also follows from taking the time derivative of the Hamil- 
tonian constraint ( p.24| ) and taking into account that S = K and T = K^, which is consistent 
with the discussion in chapter 4. 

We can actually verify that the following relations have to hold: 



Ay 



4 + 2(z/' + -)t2-e^'^^^^t, 



1 

t; - e'^v'ti + ( 3z/' + - Vs - ^ ( l{l + l)(rt6 - h) + 2t^ 



t'e + 2 ( z/' + - ^ ( (1 - /(/ + l)% + tg 



(5.28) 
(5.29) 
(5.30) 



The constraint equations serve as initial value equations that have to be solved on the initial 
slice t = in order to obtain valid initial data. Unfortunately, there is no unique way to 
solve those equations. This is due to the fact that to a particular solution of the inhomogeneous 
equations, we can always add a solution of the homogeneous equation, which would correspond 
to adding some arbitrary gravitational waves. The problem of finding the "right" initial data 
that represent only the perturbations which are due to the presence of the particle and contain 



no additional radiation will be discussed in more detail in section ^ 
For a particle falling from rest we would have v(t 

t2 



0) = and L = 0, and therefore 



= 0. Thus the momentum constraint ( |5.27D can be trivially satisfied by setting the 
extrinsic curvature variables to zero, which corresponds to time symmetric initial data. We then 
are left with solving the Hamiltonian constraint ( |5.24| ). Of course, a particle falling from rest 
would fall radially towards the neutron star and eventually hit its surface. Since we want to 
avoid such an impact, we have to give the particle some angular momentum. We also will give 
it some radial velocity in order to decrease the time it takes to come close to the neutron star. 



However, this means that we have to solve the momentum constraint (5.27), too 



5.2 Numerical implementation 

In order to solve the equations on the computer, we have to use the explicit forms of the spherical 
harmonics Y/m. The perturbation equations without particle are degenerate with respect to m, 
since the background metric is spherically symmetric. However, the presence of the particle 
breaks this symmetry and we have to consider the various m-cases. Fortunately, we do not have 
to consider all possible values of m for a given value of /, since for negative m the spherical 
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harmonics just undergo sign change and phase shift (1^^ = (— l)™!/ The advantage of 
putting the particle in the equatorial plain (6 = |) is that in the polar case we only have to deal 
with multipoles with m = I, 1 — 2, the remaining ones with m = I — 1, 1 — 3, ... are axial 
multipoles. Since the evolution code can only handle real valued perturbations, we have to treat 
the real and imaginary parts of the spherical harmonics separately. Finally, all the equations will 
be solved on a finite grid, hence we have to approximate the (5-function by a narrow Gaussian 

1 (r-fl(t))^ 

6{r — R(t)) ^ — =6 2^ , a small . 
(Tv2vr 

We also need first and second derivatives, which can be expressed as 

5'{r-R{t)) ^ ^^^^lfls{r - R{t)) 



S"{r - R{t)) ^ 1 (^-^^^^-^ -l]6{r- R{t)) . 



For this approximation to be valid, we have to ensure the convergence of the solution for a — > 0. 
This can be done in two different ways. On the one hand we can look at the convergence of 
the waveforms that are obtained in the evolution, and on the other hand we can monitor the 
violation of the constraints. A possible way do to so is to monitor the following quantity 



re 



„2u 



(rhs of (|5^)) dr , (5.31) 



where the domain of integration is the region outside the neutron star. In the limit a ^ 
it is J = 1 throughout the whole evolution, which is, of course, not true in the discretized 
form. Numerically, we cannot take this limit with a fixed grid size, since eventually we cannot 
sufficiently resolve the Gaussian. Therefore, by performing this limit, we also have to decrease 
the grid spacing in order to resolve the narrowing Gaussian. We do so by keeping the ratio a / Ax 
constant throughout the sequence. From the numerical data we find second order convergence. 
In the following, we will use a/ Ax ~ 0.15, which gives us a decent resolution of the Gaussian 
and its derivatives. 

There is still a subtle point. In deriving the source terms (see Appendix C), we have tacitly 
transformed the particle coordinate R into the spacetime coordinate r since the presence of the 
5-function makes a distinction unnecessary. However, in the evolution equations we have to 
take derivatives of the source terms with respect to r but not with respect to R, and therefore 
we would obtain different source terms if we had not changed the R's into r's. As an example 
consider the following two source terms 

S,{r) = f{r)6{r-R) (5.32) 

and 

S^ir) = fiR)6{r - R) , (5.33) 
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which are equivalent because of the presence of the (5-function. But if we now differentiate 
5*1 and 5*2 with respect to r we obtain for 5*2 just the derivative of the 5-function, whereas for 
5*1 we also have to differentiate /. Analytically this does not make a difference, but if we 
approximate the 5-function by a Gaussian, then the two expressions for 5*1 and 5*2 and their 
respective derivatives are different. To gain accuracy, we should have kept e.g. an expression 
like r'^dcj)/ dT5{r — R) as Lr"^ / B^5{r — R) and not just as L5{r — R). However, this is rather 
cumbersome, and for the numerical evolutions the actual difference is negligible, so we have 
assumed the source terms to be of the form of ^i. 



5.3 Setting up the initial conditions 

As we already know from chapter 4, any construction of initial data for a particle falling from 
rest or starting with some initial velocity involves solving the Hamiltonian constraint ( |5.24[ ). 



However, ( |5.24| ) is one equation for the two metric quantities S and T, which means that we 
have more or less the same freedom in choosing the initial values as in chapter 4. Let us denote 
the initial values of S and T at t = by So and Tq, respectively. We then can, for example, 
either choose Tq = and solve for Sq, or do it the other way round and set 5*0 = and solve for 
Tq. In the former case, we would have to solve a first order equation for 5*0, the latter case would 



lead to a second order equation for Tq. In Figs. ^]T] and we show those two possibilities. By 



comparing the different shapes of 5*0 and Tq we might already deduce that the case in which we 



set To = and solve for Sq (Fig. pTT| ) is less favorable since Sq exhibits a discontinuity at the 



location of the particle. In the other case (Fig. Tq still exhibits a kink, but it is continuous. 



To assess which choice is the more natural, we consider a particle initially at rest in flat 
spacetime. Of course, the particle will remain at rest since there is no matter around that could 
attract the particle. Thus, the perturbation of the spacetime that is created by the particle will be 
stationary. Hence the equations of motion for the metric perturbations S and T will read ^ = 
and ^ = 0. From this conditions and from the Hamiltonian constraint it then follows that S 
has to vanish. 

Of course, in the presence of the neutron star, those arguments do not hold any more, and 
5*0 = will not be the right choice of initial data, but if the particle initially is far enough away 
from the neutron star, the error in setting So = should be very small. This error actually 
corresponds to an introduction of an extra amount of gravitational radiation, which is not at all 
related to the radiation that is emitted when the particle moves through the spacetime. This extra 
amount will propagate during the evolution and eventually excite oscillations of the neutron 
star. However, if we put the particle far enough away from the neutron star, the strength of 
the induced oscillations should be small compared to the ones excited when the particle comes 
close to the neutron star, which can be confirmed by the numerical evolutions. Close to the 
neutron star, where the gravitational field is strongest, setting So = will cause S to bulge up 
and send a wave towards infinity. Again, this amount of radiation is by far smaller than the 
radiation that will come directly from the particle itself. 

In Figs. through we show the evolution of the two possible choices of initial data 



for a particle falling from rest. In Fig. 5.3 and Fig. 5.5 we show the evolution of S and T for 
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Figure 5.1: Profile of So with Tq 
location. 



0. Note that So exhibits a discontinuity at the particle's 




160 



Figure 5.2: Profile of Tq with 5*0 = 0. Note that Tq is continuous in contrast to the initial profile 



of 5*0 in Fig. 5.1 



the initial data of Fig. pTT] , and in Fig. |5]^ and Fig. ^]6| the evolution of S and T for the initial 
data of Fig. |5]^. The differences are obvious. The initial shape of T in the latter case is almost 
unchanged during the evolution, whereas S starts to aquire its "right" shape. In the other case 
we can see a huge burst of radiation propagating in both directions. In the same time T is 
acquiring its "right" shape. We also include the evolution of the Hamiltonian constraint, where 
we plot the righthand side of ( |5.24[ ), which monitors the "path" of the particle. In both cases 
the graphs agree as it should be, of course, for the Hamiltonian constraint should only monitor 
the energy density of the particle, independent of any gravitational waves. Lastly, we also show 
the metric function S and T after a certain time t, which clearly shows that regardless of the 
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chosen initial data they will adjust themselves to their proper values after having radiated away 
the superfluous initial wave content. 

We should note that we cannot escape the whole ambiguity of how to choose the variables S 
and T by using the Zerilli formalism instead. There, any regular initial data are valid initial data 
since by construction the Zerilli function always satisfies the constraint equations. Hence, we 
cannot see a priori whether ot not the chosen initial data will have additional radiation content. 

If the particle initially is not at rest, in addition to solving the Hamiltonian constraint ( |5.24D 



we also have to solve the momentum constraints (|5.25a|) and (|5.25bp. Of course, as with the 



Hamiltonian constraint we again are faced with the same kinds of ambiguity in solving that 
equation. 

Now, for a particle initially being at rest or very slow, or for circular orbits, we do not have 
to care about what kind of initial data we choose, for we do not have to worry about the initial 
gravitational wave pulses, since they travel with the speed of light and will long be gone when 
the particle, which is much slower, comes close to the star. However, if the particle's initial 
velocity is close to the speed of light, then the particle "rides" on its own wave pulse and it will 
be not clear any more whether the excitation of some particular modes of the neutron star is 
due the the particle itself or due to the initial burst, which comes from the inappropriate initial 
data. If the particle is slow enough, those two effects can clearly be distinguished. For very fast 
particles this is not possible any more and this is particularly bothersome, insofar it is known 
that a pulse of gravitational waves will predominately excite w-modes. Now, if we have a wave 
signal from a particle that grazes a neutron star with almost the speed of light, and we find that, 
indeed, there are some traces of a w-mode, how then can we make sure that this is a "real" 
signal and not an artefact due to the inappropriate initial data? 

To obtain an approximate answer, we again turn to the flat space case. In this limit, the 
equations governing the evolution of S and T reduce to two simple coupled wave equations 
with a source term that takes the presence of the particle into account: 

S + 167ifxE—5ir - R{t))Yi*^ (5.34a) 





_ d^S 


l{l + l) 












m + i) 









T + AS- %Ti-^5{r - R{t))YC^ . (5.34b) 

In flat space the particle moves on a straight line with constant velocity v, hence it is 

R{t) = ro + vt , (5.35) 

with ro being the initial location of the particle. Furthermore, the normalized energy E is just 
the Lorentz factor 

E = -^=L= . (5.36) 



V 



2 



It is interesting to note that the wave equation for S (|5.34aD is totally decoupled from the one 



for T. However, the solutions of ( |5.34aD and ( |5.34b| ) have to satisfy the flat space Hamiltonian 
constraint, which reads 



dr"^ dr V 2 / r 
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Figure 5.3: Evolution of S for To = and So from Fig. pTT] . We can see two bursts of gravita- 
tional waves that propagate both in and outwards. The ingoing pulse gets reflected at the origin 
and travels back outwards again. 
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Figure 5.4: Evolution of S for S'o = and To from Fig. Here, we can see a wave emerging 
from the vicinity of the neutron star and travelling outwards. However, the amplitude is much 



smaller than in Fig. 5.3 
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Figure 5.5: Evolution of T for Tq = and So from Fig. pTT] . We can see the build-up of T, 
which is disturbed by the reflected part of the wave that was sent out by the particle. 
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Figure 5.6: Evolution of T for So = and Tq from Fig. |5^ . It is evident that T basically has 
its "right" shape right from the beginning. 
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Figure 5.7: Evaluation of the Hamiltonian constraint during the evolution for Tq = 0. The 
Gaussian shape of the particle is clearly visible and has been chosen to be very broad for a 
better visualization. The particle initially is at rest and starts falling towards the neutron star. 
On the very right side we can see the matter oscillations of the neutron star. 
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Figure 5.8: Evaluation of the Hamiltonian constraint during the evolution for Sq = 0. As 
expected, it basically does not differ from Fig. ^77[ for the gravitational waves do not give any 
contribution. Note, however, that the matter oscillations of the neutron star are much smaller 
than in Fig. 



5.3. Setting up the initial conditions 



101 




20 40 60 80 100 



r in km 

Figure 5.9: Plot of the variables S at the end of the evolution of the different initial data. The 
artificial radiation of the initial data has been radiated away and S has assumed its "right" profile 
that is independent of the initial data. The difference between the different profiles is due to the 
fact that the initial data with Tq = contain much more radiation, which excites the neutron 
stars to pulsations, which in turn disturb the profile of 5*. 
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Figure 5.10: Plot of the variables T at the end of the evolution of the different initial data. 
Here, too, the artificial radiation of the initial data has been radiated away and T has assumed 
its "right" profile, which is independent of the initial data. In contrast to S there is almost no 
difference in the two profiles. 



We now seek for an exact solution of ( |5.34a| ) that obeys the right boundary conditions at the ori- 
gin and at infinity. Once it is found, we may use ( ]5.37D to numerically compute the appropriate 
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T. We state that a solution for ( |5.34aj ) is given by the following series ansatz: 

/ °° ^.2i+l+l 



\l-2i-2 



(5.38) 



i=0 ^ 



where 6 is the Heaviside function, which satisfies 

^ / X 1 0, X < 

Q(x) = < ' . (5.39) 

\l, a;>0 

Continuity at r = + vt requires that 

oo oo 

= (5.40) 

i=0 i=0 

The overall amplitude will be determined by A, hence we deliberately may set 

oo oo 

J2a, = = 1- (5.41) 

i=0 i=0 

The amplitude A and the coefficients and bi can be found by plugging ( |5.38| ) into ( |5.34a| ). 
For A we find 

A = ^6vr^E3,.^^ 

2/ + l + 2E.^o2(a, -6,) ' 
and the coefficients and bi are determined by the following recursion relations 

, (2^ + / + 3)(22 + / + 4) 
2{i + l)(2z + 2/ + 3) 
^ ^il-2t-2){l-2t-3) ^^^^^ 
= 2(. + l)(2.-2/ + l) • ^'-''^ 

It is interesting to note that whereas the series in will extend to infinity, the series in bi will 
always terminate because one of the two factors in the numerator will become zero for some i. 
The series in will converge if and only if |f | < 1. 

In Fig. |5.1 1| we show the initial data obtained from ( |5.38| ) for a particle that is located at r = 
500 km with different initial velocities. Here, we plot rS and T/r since only those quantities 
can be meaningfully compared with each other. For t> = it is rS* = but the amplitude of rS 
grows rapidly when the particle's velocity approaches the speed of light, whereas the peak of 
T/r slightly decreases. In the ultra-relativistic limit rS totally dominates over T/r. 

Of course, the above solutions for S and T are not valid any more if we consider the parti- 
cle in curved spacetime, because they would violate the constraints. However, we can still use 
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( |5.38p as a prescription for S and then use the curved space Hamiltonian constraint ( |5.24| ) to 
solve for T. Furthermore, we can compute K from K = dS/dt and then use the momentum 
constraint ( |5.27D to compute K^. As long as the particle does not have any angular momentum 
and is far away from the neutron star, the thus obtained initial values should be a good approx- 
imation for a boosted particle on a Schwarzschild background. However, if the particle has a 
large angular momentum, there will be additional source terms in the evolution equations and 
our approximation should break down. However, we are mainly interested in trajectories which 
come very close to the neutron star, and hence have only small angular momentum. In this case 
the above prescription still yields good initial data. 



5.4 Numerical results 

To compare the code with known results, we consider a particle in a circular orbit with radius 
ro around the neutron star and compute the radiated energy at infinity. Numerically, this can be 
accomplished by evaluating the Zerilli function Z, which can be computed from S and T by 
means of formula ( |4.33| ), at some large distances. 

The radiated power for some particular values of / and m can then be computed from 

dEi^ _ 1 ([+2)! . 2 ..... 

(For a derivation, see Appendix C). Cutler et al. have numerically computed the normalized 
gravitational power {M/ nYdEim/ dt that is radiated by a particle orbiting a black hole. In table 
II they compile the multipole components for tq/M = 10. To compare the output of our code 
to those results, we choose the mass of the neutron star to be M = 1.99 km and the orbit of 
the particle to be at tq = 19.9 km in order to obtain the ratio of r^/M = 10. The mass of 
the particle is set to /i = 1 km. In Fig. (|5.12[) we show Zim as a function of time extracted 



at r = 500 km for / = m = 2. After some wave burst that comes from the inappropriate 
choice of initial data, we see that the signal is periodic with a frequency of twice the orbital 
frequency of the particle. The amplitude is about 0.0076, which corresponds to a radiated 
power of (M/ iiYdE22/dt = 5.46-10"^, which is in excellent agreement with Cutler et al. ,who 
obtain a value of (M/ iif'dE22/dt = 5.388-10^^. The slightly higher value of our result may be 
due to the fact that at r = 500 km the Zerilli function is still somewhat off its asymptotic value, 
which will be a little bit smaller. Also, a black hole absorbs some of the radiation, whereas a 
neutron star will re-emit all of the incoming radiation. We find that we agree with all the polar 
modes compiled in table II of Cutler et al. within a few percent, only for the case / = 3 and 
m = 1 we find the radiated power to be 5.9- 10^^° instead of their value of 5.71-10"^. 

It is clear that a particle on a circular orbit will not excite the eigenmodes of the neutron 
star in a significant manner, unless the orbital frequency is close to the frequency of a stellar 
quasi-normal mode. In this case, we can have a resonant excitation of this particular mode 
and the radiated energy flux can drastically increase []62[]. It is clear, however, that only the 
low-frequency modes like the /-mode can be excited by this mechanism; the frequencies of the 
w-modes are much too high to lie in the frequency range of the orbiting particle. 
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Therefore the only way to possibly excite w-modes are very eccentric orbits, where the 
periastron is very close to the surface of the star, or scattering processes with very small impact 
parameters. The investigation of the latter is the main objective of this work. 

Before presenting the results for scattering orbits, we should comment on some difficulties 
that are related with the use of our evolution code. 

As we have learned in the previous chapter, the computation of the Zerilli function is some- 
what troublesome, for the violation of the Hamiltonian constraint during the numerical evolu- 
tion prevents the exact cancellation of the growth in T. However, we managed to go around 
this problem by extracting the Zerilli function Z very close to the star and then using the Zerilli 
equation itself to propagate Z towards infinity. This lead to very reliable results. 

There should be nothing that prevents us from doing the same thing for the particle. Of 
course, we would have to include the appropriate source term which takes care of the presence 
of the particle. Unfortunately, the inclusion of this term is not straightforward, but for our 
purposes it is not really necessary either. 

Our concern is not so much to obtain the quantitative amount of energy that gets radiated 
away in a scattering process. We are much more interested in some qualitative statements 
concerning the relative excitations of the various neutron star modes, and especially whether or 
not the particle can excite the w -modes. For this purpose it is enough to look at the waveforms 
of S and T, which also contain all the relevant information. 

Of course, it is not totally impossible to compute a quite accurate Zerilli function as we 
have already demonstrated for the circular orbits. Here, the problems with the amplification 
of the high frequency components are not present, since the orbital frequency of the particle is 
comparably low. Still, we were forced to resort to quite high resolution in order to obtain the 
desired accuracy. 

Before we go on with the discussion we should clarify one point that might give rise to 
some confusion. We have decomposed the particle, which is represented by a 3-dimensional 
5-function, into its various multipole components, which are labeled by / and m. Each of those 
multipoles represents an infinitely thin matter shell, whose surface density is proportional to the 
spherical harmonic Yim(^, 0). It is only by the summation over all the shells, i.e. over all I and 
m, that the particle gets localized at the particular point (R, 0, $) in space. In the following we 
focus on / = m = 2, but we will still use the term "particle" though we better had to speak of a 
"quadrupole shell". 

Ideally, we should extract the wave form at r = cxd, which is clearly not possible on a finite 
grid. We have to record the signal at some finite distance from the star. But if the particle moves 
on an unbounded orbit, it means that it will eventually cross the location of the observer. 

Now, the particle is always slower than the propagation speed of the gravitational waves, 
which, of course, propagate with the speed of light. Hence, the observer will usually see the 
wave signal before he is hit by the particle (better: quadrupole shell). However, when the 
particle is very fast and the observer is not far enough away from the neutron star, not enough 
time has elapsed for the interesting part of the signal to cross the observer before the particle 
arrives. That means that the signal will be a superposition of the "real" gravitational wave signal 
and the influence of the gravitational field of the particle itself. 

The farther the observer moves outwards, the better the separation of the two components 
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in the signal can be made. Furthermore, the influence of the particle, i.e. the strengths of the 
source terms in the equations for S and T decrease with for S and with 1/r for T (see 
equations ( |5.34D ), whereas the amplitude of the outgoing signal remains constant for S and 
grows linearly for T. But this means that the presence of the particle affects the wave form of 
T much stronger than that of S. 

Of course, the effective interfering influence of the particle strongly depends on the actual 
excitation strength of the neutron star oscillations. The latter in turn depends on the value of the 
particle's turning point rt, which is the smallest distance between the particle and the neutron 
star. The closer the particle approaches the neutron star, i.e. the smaller rj is, the higher is the 
amplitude of the induced oscillations. For large the induced fluid oscillations are so weak 
that they will totally be buried within the gravitational field of the particle. This is depicted 
in Fig. |5.13| , where we show the wave forms of S and T for two different orbits with turning 
points of rt = 10 km and = 50 km. The influence of the gravitational field on the wave forms, 
especially on T, is obvious. 

In Figs. [5.14| through |5.16| we show the wave forms of S and T for different initial radial 
velocities vq and different parameters of the turning point rt. Figures |5.17| through |5.19| show 
the corresponding Fourier transformations for the orbit with rt = 10 km. 

From Figs. |5.14| and |5.15| and their corresponding power spectra Figs. |5.17| and |5.18| it is 
clear that for "small" initial velocities vq < 0.5 there is no hint of a u;-mode excitation in the 
signal. Instead the signal consists of a first pulse, which comes from the part of the particle's 
orbit close to the turning point, where the particle is closest to the neutron star and radiates the 
most strongly. The power spectrum of this part of the signal should peak at about twice the 
angular velocity cu = d^/dt of the particle at the turning point rt, since the particle's source 
term is proportional to cos{mujt). For vq = 0.5 and rt = 10km the peak is around 5.5kHz. 

As soon as the particle leaves the star, its radiation strongly decreases and the fluid modes 
of the excited neutron star take over. Here, we find that almost all the energy is in the /-mode, 
only a tiny fraction is in the first p-mode. 

It is only for very high initial values of vq and very small values of rt that there is indeed 
a significant excitation of w-modes. By sampling different initial velocities we find that for 
Vq ~ 0.7 we can spot a trace of the first w-mode. And for vq approaching 1, we can obtain a 
quite strong excitation of the first w-mode, which can be seen in Figs. |5.16| and|5.19. 



5.5 Discussion 

With the particle we have a physical even though probably a quite unrealistic process, which 
induces oscillations in a neutron star. By unrealistic we mean that this mechanism is unlikely to 
generate gravitational waves that are strong enough to be detectable on Earth. 

Whereas all the considered kinds of initial data of Chapter 4 could excite the /-mode and 
some of the first p-modes, this was not the case for the if-modes. There, only a special class of 
initial data was able to excite those modes. 

Andersson and Kokkotas have shown how to extract the physical parameters of a neu- 
tron star, like mass M and radius R and possibly the equation of state, if the complex frequencies 
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of the /-mode and the first w-mode are known. Of course, it is only possible to detect the first 
w-mode of a neutron star if it gets excited by some physical mechanism. 

We used an orbiting mass yU in the particle limit yU <^ M to assess whether or not it is 
possible to excite the first w-mode. We found that in order to excite a significant amount of 
if-modes the particle's velocity at infinity must be much more than 70 percent of the speed of 
light. Moreover, the excitation strength of the modes rapidly decreases as the particle's turning 
point rt increases. This is in agreement with the results of Ferrari et al. [^, who find that no 
significant w-mode excitation is observable. However, they do not consider the cases where the 
particle has an energy E ^ 1, which is necessary if it were to excite any w-modes. 

The question now arises: Can we infer anything about what happens in an astrophysical 
event? Our results show that the particle initially has to be incredibly fast in order to excite a 
significant amount of u'-mode content. It is clear that there do not exist any astrophysical events 
at all that could accelerate the particle (which, of course does not represent an elementary par- 
ticle, but some heavy extended object, like a planet or an asteroid or even another neutron star) 
up to more than half the speed of light. Hence, we might safely conclude that any astrophys- 
ical scenario that can be simulated by a particle orbiting a neutron star will not produce any 
detectable amount of w-modes. 

But what if the object collides with the neutron star? It has been suggested [ |65l ] that in 
this case there might be some significant excitation of ty-modes. In particular, simulations of a 
particle that radially falls onto a neutron star were done by Borelli [ |64| ] who found that, indeed, 
both fluid and w-modes were excited to a significant level. However, one serious drawback 
of those calculations is that as soon as the particle hits the surface of the star, the calculation 
is stopped. This leads to an overestimation of the high frequency components in the power 
spectrum. Besides, in this study only axial modes were investigated, it is therefore clear that the 
w-modes will always show up in the spectrum because they will not be screened by the presence 
of fluid modes, which can happen in the even parity case. Now, even in the particle limit, it is 
not possible to simulate a collision since it not clear at all what happens when the particle hits 
the surface. A first step might be, though, to let the particle go right through the neutron star 
without being affected by the presence of the neutron star matter. In case there is a significant 
amount of ly-modes one might conclude that in a realistic scenario, where the physics of the 
impact is included, there still might be w-modes present, maybe they are even more strongly 
excited than in the unphysical "tunnelling" case. However, our view is that even in the impact 
case the particle still must have an initial velocity that is a significant fraction of the speed of 
light. Since for scattering orbits with initial speeds less than 50 percent of the speed of light 
there is no hint of any w-mode presence even for the closest trajectories, we believe that this 
should not change very much even if the particle's turning point lies inside the neutron star. Of 
course, in this case it seems reasonable to infer that the p-modes of the fluid would be much 
more strongly excited. 

Our results in chapter 4 indicate that the presence of the if -modes in the signal is closely 
related with the value of S. This particular role of S seems also to be confirmed in our par- 
ticle simulations, for we have seen that for ultra-relativistic particles the value of rS greatly 
dominates the value of T/r. And it is only there that we find a significant amount of w-modes. 

Still, the question remains as to what happens in a realistic astrophysical scenario. We 
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have excluded the possibility of w-modes excitation by means of a realistic scattering process. 
However, it is not possible to confer our results to the merger process of a binary neutron star 
system, for the latter cannot be adequately described within perturbation theory. Only if the 
final object does not immediately collapse to a black hole, it might wildly oscillate and emit 
some significant radiation through w -modes. 
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Figure 5.1 1: Initial data for a particle with different initial radial velocities v, where S is given 
by ( |5.38[ ) and T is obtained by solving ( |5.37| ). For small velocities, T/r dominates rS, whereas 
for ultra-relativistic velocities, rS dominates T/r. 
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Figure 5.12: Evolution of the real and imaginary parts of Z22 at r = 500 km. After the initial 
radiation bursts the wave forms show periodical oscillations with twice the orbital frequency. 
The amplitude is about 0.0076, which corresponds to a radiation power of {M/ nydE22/dt = 
5.46-10-^ 
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Figure 5.13: It is shown tlie difference of excitation strengths for the two different grazing radii 
rt = 10 km (solid line) and = 50 km (dashed line). The waveform of S (upper panel) for 
rt = 10 km is totally unaffected by the presence of the particle, whereas for = 50 km the 
gravitational field shifts the amplitude of the signal to higher values. For T (lower panel) this 
effect is much more pronounced and can already be detected for rt = 10 km. For rt — 50 km, 
the oscillations of the neutron star are totally buried in the gravitational field of the particle. At 
about t — 10 ms the particle crosses the observer who is located at robs — 1000 km. 




Figure 5.14: Wave forms of S and T for the three different grazing radii rt = 
and rt — 30km. The initial velocity of the particle is i^o = 0.1. 
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Figure 5.15: Wave forms of S and T for the three different grazing radii = 10km, = 20km 
and Vt = 30km. The initial velocity of the particle is = 0.5. 
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Figure 5.16: Wave forms of S and T for the three different grazing radii Vt = 10km, rj = 20km 
and rt — 50km. The initial velocity of the particle is vq — 0.97. 
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Figure 5.18: Fourier transformation of the wave forms for = 0.5. 
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Figure 5.19: Fourier transformation of the wave form for = 10km and Vq — 0.97. The scales 
are arbitrary. The presence of the first w-mode is clearly visible. 



Appendix A 



Relations between the Regge-Wheeler and 
Mathews-Zerilli harmonics 



The evolution equations will be expanded using the Regge-Wheeler harmonics [3^/^]^,^, 
which are defined in [|T7p, and which do not form an orthonormal set. However, in order to 
perform the expansion of the energy-momentum tensor (|5.1|), we need an orthonormal set. This 
is given by the Mathews-Zerilli harmonics [3^/^]^^ [ [74| , |23| ], which are orthonormal with respect 
to the following inner product 
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[yi'm']lJ-i^^^ = SaSaA'Sii'S„ 



where the asterisk denotes complex conjugation and 



[ytr 



(A.l) 



(A.2) 



Because of the use of the inverse Minkowski metric r/^'^ to raise the indices, the inner product 
dA.lD is not positive definite and it is 
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For the sake of completeness we list the whole set of the 10 tensor harmonics 3^f^, where we 
correct some errors in [ |6T| ] : 
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and Yim are the ordinary scalar harmonics. There are three odd parity harmonics y^^, yj^ and 
yl^, the remaining seven ones have even parity. 

In terms of the Mathews-Zerilli harmonics y^ the Regge-Wheeler harmonics y(^ read: 



yt = -i\/25?L 

yim yim 

yfm = l^/W^yl 

yL = ly^i(7TT)5^,; 
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The inverse relations are given by 
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Appendix B 

The source terms of the particle 



The particle's energy-momentum tensor ( p7T| ) will be expanded using the orthonormal Mathews- 
Zerilli harmonics defined in appendix A: 
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1=0 m=-l A=l 



By using the orthonormality condition ( |A.1[ ) we can compute the coefficients through 



tA = / [ytr'"T,.dQ. (B.2) 
We thus obtain the following set: 
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Here, F^^, IV^^ and X^^ are functions of the particle's angular position 6 and $ parametized by 
the coordinate time t, therefore all derivatives with respect to proper time r are to be understood 

as 
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Since for the decomposition of the evolution equations we use the Regge-Wheeler harmonics 
y^, which do not form an orthonormal set, we have to expand the energy-momentum tensor 



( p.l[ ) in Regge-Wheeler harmonics as well: 
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Of course, we cannot obtain the Regge-Wheeler coefficients by means of a formula sim- 
ilar to ( p.2[ ), but we can construct them from the Mathews -Zerilli coefficients t^. Using the 
relationship between the different sets of tensor harmonics given in Appendix A, we find: 



Am ±Li 



Im 



Im 
2 

Im 
3 

Im 
4 

Im 

5 

±lm 

''6 

dm 



±lm 



T2 



ir 



1.1m 



llm 



V2/(/ + l) 



ir 



vw 



118 



Appendix B. The source terms of the particle 



2r^ 

^Im £im 

' ~ v/2/(/ + l)(/-l)(/ + 2) ' 



Im _ ( flm I / /(/ + 1) 



Am llm 

'^lO ~ / ^,,, , , ^X ^IO 



v/2/(Z + l)(/ -!)(/ + 2) 
We now restrict the motion of the particle to the equatorial plane = |. In this case it is 



and sin 6 = 1, and we can use the geodesic equations ( |5.3aD and ( |5.3c[ ) to substitute all 



dB 
dr 

expressions containing derivatives with respect to proper time r. Furthermore, we can obtain 
quite simple relations for the derivatives of Yi*^: 

-^Y;^ = -imF4 (B.5) 
^Y.l = (m^ - /(/ + 1)) Y,l . (B.6) 

This then gives us a somewhat simpler set of coefficients: 
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where v{t) = ^ is the radial velocity of the particle. The field equations also require the 
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Using the explicit forms of the coefficients, we obtain 

t'- = ^6{r - Rit)) (Ev'{t)e'' -E + e'^^ Y^^ , (B.9) 
and by making use of the geodesic equation ( |5.3b| ) we can reduce this expression to 

t'" = -^^"^ - ■ (B.IO) 



Appendix C 

Derivation of the radiation extraction 
formula 



The radiated energy emitted per unit time and unit angle is given by 
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The metric components have to be in the radiation gauge, which asserts that only 
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whereas all other components must fall off faster than 0{l/r). Furthermore, it must hold that 
h[e][e] = ~h[(f)][(f)] + 0{l/r'^). We now expand the angular parts of the metric into the orthonormal 
Mathews-Zerilli tensor harmonics given in Appendix A. Since we focus on the polar perturba- 
tions, we only need to consider 3^g™ and 3^9™, and we have 



h 



E 

Lm 



1 



Afsin 



1 G^Zm^Zm 



(C.8) 
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where we use 

N:= v/2Z(/ + l)(/-l)(/ + 2) 
This gives us 



(C.9) 



dtdfl ~ 16^ 



Lm 



Lm 



Ar2 



(C.IO) 



where we have to take the square moduli because we are dealing with complex metric coeffi- 
cients. To obtain the total energy flux, we have to integrate over a 2-sphere 



dE 



r 

IQtt 



1 



l,m,l' ,m' 
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(C.ll) 



Now, because of the orthonormality relation the integral is just 



I s'2 V sin^ 9 
and we obtain 

dE 



+ Wi^wr,^, dn 



W'Jmm' ; 



(C.12) 



dt 3271 



\Glm\ 



(C.13) 



l,m 



Since the actual metric perturbations are derived using the expansion into the Regge- Wheeler 
harmonics, we have to express ( |C.13| ) in terms of those variables. The relation between Gim 
and the Regge- Wheeler variable Gim is easily derived to be 



G 



N 



lm 



-G, 



lm 



(C.14) 



which leads to the following radiation formula 



dE 
'dt 



r 

64^ 



5^/(/ + l)(/-l)(/ + 2)|G 

l,m 



lm 



(CAS) 



Note, however, that Gim is still in the radiation gauge and not in the Regge-Wheeler gauge. In 
the Regge-Wheeler gauge we have Gim = 0, and the above formula would not make any sense. 
We now have to devise a way which relates Gim to the variables in the Regge-Wheeler gauge. 
This can be accomplished as follows. Following Moncrief []2^], we construct some gauge in- 
variant quantities, in particular the Zerilli function, and examine their asymptotic behavior. We 
then will show that the Zerilli function has the same asymptotic behavior as Gim in the radiation 
gauge. Hence, the radiated power is proportional to the time derivative of the Zerilli function, 
which can be easily constructed from the metric variables in the Regge-Wheeler gauge. 
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Appendix C. Derivation of the radiation extraction formula 



Let us write the general expansion of the even parity perturbations of the spacial part of the 

metric in terms of Regge-Wheeler harmonics (Here we use the notation of Moncrief and focus 
on fixed values I and m:) 

Kr = e^^H^Yi^ (C.16) 

hre^hi^^, hr4, = (C.17) 

hee = (^KYi^ + G^p-^ (CAS) 

2 

he4, = jGXi^ (C.19) 

V = sin^ e (^KYim - G (^1(1 + l)Yi^ + ^ . (C.20) 

We can then construct the following gauge invariant quantity 

gi := 4re-^^A;2 + rl{l + l)ki , (C.21) 

where 

ki^K + e-^^ (rG' -'^hij (C.22) 



and 



h = l [e^'H, - {re'K)') . (C.23) 



We now have to determine the asymptotic behavior of qi. To do so, we have to know the 
asymptotic behavior of the metric coefficients hi, H2, K and G, which, of course, will depend 
on the chosen gauge. However, because of its gauge invariance, qi will always have the same 
asymptotic behavior, regardless in what gauge the metric variables are. 

Therefore we assume our metric variables to be in the radiation gauge, where we know the 
asymptotic behavior. Since h[o][e] — — ^[<^][(/>] + 0{l/r'^), the sum h[o][o] + h[^][4>] has to be of 
order 0(l/r^), i.e. 

horn + /^MM = {2K - 1(1 + 1)G) Yi^ = 0{l/r^) . (C.24) 
But this can only hold if 

2X = /(/ + l)G + C(l/r2) (C.25) 
for the leading order. Using this relation we find for k-i and k2 

ki ^ K + rG' -l{l + l)G + rG' (C.26) 
2 

k2 ~ -^K-^'-K' ~ -iz(Z + l)(G + rG') , (C.27) 
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where "a ~ b" means a equals b plus higher order terms. This leaves us with 

qi ~ ^/(/ + l)(/(/ + l)-2)G'. (C.28) 
Moncrief finally defines the quantity 



/(/ + !)_ 2 + Ml 



which satisfies the Zerilli equation. The leading order for Q is then 

Q ~ U{l + 1)G, (C.30) 
which is equivalent to 

This relation only holds if G is given in the radiation gauge. But in the radiation formula ( |C.15D 
this is the case and we may substitute G by means of ( |(J.31| ), which yields 

dE _ 1 ^ (/-!)(/ + 2) . 

l,m 

We now have to express Q in terms of our actual perturbation variables S and T, which are in 
the Regge-Wheeler gauge. We first recall that in the Regge-Wheeler gauge, it is 

hi = (C.33) 

H2 = -+rS (C.34) 
r 

K = - (C.35) 

r 

G = . (C.36) 
Then we compute ki and k2, which is readily done: 

ki = - (C.37) 

r 

fc. = ^e-(^ + r^)-eMe^T)') . (C.38) 
Next is qi 

q, = 2re-'^ (^e^^ + rS^ - (e^T)'^ + /(/ + 1)T (C.39) 

= 2e-^^ ((1 - rX') T + r^S- T) + /(/ + 1)T . (C.40) 
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Appendix C. Derivation of the radiation extraction formula 



Using 



A' 



M 



r (2M - r) 
and collecting terms yields 

2M 



gi = - 1 



, 2M-r(2 + /(/ + 1)) . 
2rT' + ^ \. '' T - 2r^S 



If we furthermore define Z to be 



r-2M 



/(/ + !)- /(/ + !)/(/ + !)_ 2+ ' 
we find from ( 1C3T] ) that 

G ~ — , 



and the radiation formula ( |C.32D finally reads 



l.m 



with 



2M 



/(/ + l)/(/ + l)-2 + ^ 



647r^ (^-2)! 

l,m, ^ ' 



, 2M-r 2+ +1 ^ . ' 
r-2M 



(C.41) 



(C.42) 



(C.43) 



(C.44) 



(C.45) 



(C.46) 



which agrees with our definition of the Zerilli function in ( |4.32D 
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Dynamical instability of gaseous masses approaching the Schwarzschild limit in gen- 
eral relativity. — S. Chandrasekhar, Phys. Rev. Lett. 12, 114 (1964) 

The dynamical instability of gaseous masses approaching the Schwarzschild limit in 
general relativity. — S. Chandrasekhar, Astrophys. J. 140, 417 (1964) 

Abstract: In this paper the theory of the infinitesimal, baryon-number conserving, adiabatic, radial oscillations of a 
gas sphere is developed in the framework of general relativity. A variational base for determining the characteristic 
frequencies of oscillations is established. It provides a convenient method for obtaining sufficient conditions for the 
occurence of dynamical instability. The principal result of the analysis is the demonstration that the Newtonian 
lower limit 4/3, for the ratio of the specific heats 7, for insuring dynamical stability is increased by effects arising 
from general relativity; indeed, is increased to an extent that, so long as 7 is finite, dynamical instability will 
intervene before a mass contracts to the limiting radius (> 2.2bGM / c^) compatible with hydrostatic equilibrium. 
Moreover, if 7 should exceed 4/3 only by a small amount, then dynamical instability will occur if the mass should 
contract to the radius 

K 2GM 

Rc = TTT^ (7-4/3) 

7 — 4/3 

where K is a constant depending, principally, on the density distribution in the configuration. The value of the 
constant K is explicetly evaluated for the homogeneous sphere of constant energy density and the polytropes of 
indices n = 1, 2, and 3. 

A catalogue of methods for studying the normal modes of radial pulsation of general- 
relativistic stellar models — James M. Bardeen, Kip. S. Thorne and David W. Meltzer, 
Astrophys. J. 145, 508 (1966) 

Abstract: A catalogue is presented of various methods now available for studying the normal modes of radial pulsa- 
tion of stars in general relativity. The methods are classified according to the normal-mode properties which they 
investigate; and the relative merits of the methods are discussed. 

Normal modes of radial pulsation of stars at the end point of thermonuclear evolution 

— David W. Meltzer and Kip. S. Thorne, Astrophys. J. 145, 514 (1966) 

Abstract: The lowest three normal modes of radial pulsation are studied for stellar configurations at the end point of 
thermonuclear evolution. Two alternative equations of state are assumed - the Harrison- Wheeler equation of state, 
and the Skyrme-Cameron-Saakyan equation of state. The adiabatic index is discussed for each equation of state, 
taking into account the extend to which period changes in nuclear composition are produced by stellar pulsation. 
The periods and e-folding times of the lowest three normal modes of a wide range of configurations are calculated 
to an accuracy of ±1 per cent. The eigenfunctions of some of the normal modes and the pulsation energy which 
they can store are also calculated; and a discussion is presented of the damping of neutron-star pulsations and of 
the consequent heating of the star. [...] 



Radial oscillations of zero-temperature white dwarfs and neutron stars below nuclear 
densities — G. Chanmugam, Astrophys. J. 217, 799 (1977) 

Abstract: The radial oscillations of zero-temperature degenerate stars with subnuclear densities have been calculated 
for the fundamental and first two excited modes using recent equations of state. The calculations were made under 
the assumption that (a) the dynamical time scale r 3> tjv, where rjv is the nuclear relaxation time, in which case 
the adiabatic index 7 = ■je (calculated along the equation of state) and (b) t <g; rjv with 7 = ■yc calculated at 
constant composition. For case a, it is shown that white dwarfs are dynamically stable for pc < 1.2 X lO^g cm~^ 
while neutron stars are dynamically stable for pc > 10^* g cm~^. The instability for white dwarfs in case a arises 
from principally from the phase transitions than from the effects of general relativity, confirming the conjecture of 
Baym et al. Low-density neutron stars have for the fundamental (unstable) mode "e-folding" times with a factor 
100 less than the typically expected dynamical time scale. It is shown that this is consistent with the unusual 
decay of the radial eigenfunction, a qualitative explanation of which is given. For denstities below the neutron drip 
point, it is shown that 7c = 7b except at a finite number of phase transition points. For densities just above the 
neutron drip point, it is shown analytically that -yc must be greater that 4/3 in the critical region where 7^ drops 
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below 1/3. For case 6, white dwarfs are dynamically stable for pc < 4 X 10^" g cm while neutron stars are 
dynamically stable for pc > 7 X 10^^ g cm~^ as had been suggested previously by the author. Case b is shown to 
be more appropriate for infinitesimal perturbations. It is also pointed out that arguments given earlier to indicate 
that high-density white dwarfs and and low-density neutron stars may not exist, because of instabilities over long 
time scales, are lacking rigor. 



The radial oscillations of neutron stars — Edward N. Glass and Lee Lindblom, Astrophys. 
J. 53, 93 (1983) -|- Erratum, Astrophys. J. Suppl. 71, 173 (1989) 

Abstract: We continue our investigation of the influence of the equation of state of matter (in the nuclear density 
regime) on the dynamical behavior of neutron stars. The properties of the two lowest frequency radial oscilla- 
tion modes of neutron stars, constructed from 13 equations of state, are presented. We find that the frequencies of 
both radial oscillation mods studied here are generally higher than the frequency of the corresponding quadrupole / 
mode. 



Radial oscillations of neutron stars and strange stars — H.M. Vath and G. Chanmugam, 
Astron. Astrophys. 260, 250 (1992) 

Abstract: We have determined the two lowest frequency radial oscillation modes constructed for several equations 
of state. The results for neutron stars differ from those of Glass & Lindblom (1983). We find that the central 
density where uJq = agrees in our calculations with that where the star reaches its maximal mass, as expected, but 
contrary to their calculations. The oscillation frequencies of strange stars have qualitatively a different dependence 
on the central density compared with the case of neutron stars, as the periods of all modes go the zero when the 
central density of the strange star approaches its smallest possible value. 



Radial modes of rotating neutron stars in the Chandrasekhar- Friedman formalism — 
B. Datta, S.S. Hasan, P.K. Sahu, A.R. Prasanna, I.J.M.P. D 7, 49 (1998) 

Abstract: Eigenfrequencies of radial pulsations of "slowly" rotating neutron stars are calculated in a general rela- 
tivistic formalism given by Chandrasekhar & Friedman. It is found that the square of the frequencies are always a 
decreasing function of the cental density of the neutron star. The decrease of the squared frequency is sensitive to 
the equation of state of neutron star matter, and is illustrated using realistic models. 



Radial pulsations and stability of protoneutron stars — D. Gondek, P. Haensel, and J.L. 
Zdunik, Astron. Astrophys. 325, 217 (1998) 

Abstract: Radial pulsations of newborn neutron stars (protoneutron stars) are studied for a range of internal tem- 
peratures and entropies per baryon predicted by the existing numerical simulations. Protoneutron star models 
are constructed using a realistic equation of state of hot dense matter, and under various assumptions concerning 
stellar interior (large trapped lepton number, zero trapped lepton number, isentropic, isothermal). Under prevailing 
conditions, linear oscillations of a protoneutron star can be treated as adiabatic, and evolutionary effects can be 
neglected on dynamic timcscalc. The dynamic behavior is governed by the adiabatic index, which in turn depends 
on the physical state of the stellar interior. The eigenfrequencies of the lowest radial modes of linear, adiabatic 
pulsations are calculated. Stability of protoneutron stars with respect to small radial perturbations is studied, and 
the validity of the static criteria is discussed. 



Papers on Dipole Oscillations of Neutron Stars 



Non-Radial pulsation of general-relativistic stellar models. V. Analytic analysis for 

1 = 1 — Alfonso Campolattaro, Kip S. Thorne, Astrophys. J. 159, 847 (1970) 

Abstract: The theory of small, adiabatic, dipole perturbations of a star away from hydrostatic equilibrium is de- 
veloped within the framework of general relativity. The analysis is linearized in the perturbation amplitudes. The 
"odd-parity" perturbations describe differential rotation, while the "even-parity" perturbations describe pulsation. 
In the pulsation case there are two degrees of freedom associated with the fluid motion and no degrees of freedom 
in the gravitational field (no dipole gravitational waves). During the pulsation the star's external gravitational field 
is completely unperturbed, to first order in the amplitude, despite the fact that its surface performs a finite back- 
and-forth motion. In both the pulsation and the rotational cases, the Einstein field equations are put into simple 



2 



forms suitable for numerical integration. 



A variational principle and a stability criterion for the dipole modes of pulsations of 
stellar models in general relativity — Steven L. Detweiler, Astrophys. J. 201, 440 (1975) 

Abstract: A variational principle is found for the frequencies of the dipole modes of pulsations of spherical stellar 
models in general relativity. This principle reveals that if the relativistic Schwaj-zschild discriminant, S, is not 
positive definite, then a stellar model is unstable to dipole perturbations. In addition, if S is positive definite the 
eigenfrequencies of all of the normal modes are real. 



The dipole oscillations of general relativistic neutron stars — Lee Lindblom, Randall J. 
Splinter, Astrophys. J. 345, 925 (1989) 

Abstract: We investigate the theory of the dipole oscillations of fully relativistic neutron stars. The equations gov- 
erning these modes are reduced to a third-order system; and the variational principle for the frequencies of these 
modes is reduced to a form involving only two independent functions. The equations that determine the damping 
of these modes, which is due to viscous dissipation, are presented. The frequencies and viscous damping times for a 
range of realistic neutron star models based on a number of equations of state and a range of masses are determined 
by solving these dipole oscillations equationsnumerically. 



The accuracy of the relativistic Cowling approximation — Lee Lindblom, Randall J. 

Splinter, Astrophys. J. 348, 198 (1990) 

Abstract: We investigate the accuracy of two versions of the relativistic Cowling approximation for the oscilla- 
tions of general relativistic stellar models: one by McDermott, Van Horn and SchoU, and a refinement by Finn. 
The oscillation frequencies and the viscous damping time scales of several of the lowest frequency dipole p-modes 
in the approximations are compared to the values for the quantities. We find, as expected, that the accuracy of the 
Cowling approximation improves as the number of nodes in the cigcnfunctions of the mode increases. Furthermore, 
the McDermott, Van Horn and SchoU version of the Cowling approximation is more accurate than Finn's refined 
version for the low-order p-modes studied here. 



Papers on Nonradial Oscillations of Neutron Stars 



Non-radial pulsation of general-relativistic stellar models. I. Analytic analysis for I >2 

— Kip S. Thornc, Alfonso Campolattaro, Astrophys. J. 149, 591 (1967) -|- Erratum Astrophys. 
J. 152, 673 (1968). 

Abstract: The theory of small, adiabatic, non-radial perturbations of a star away from hydrostatic equilibrium is 
developed within the framework of general relativity. The unperturbed equilibrium configuration is an arbitrary, 
non-rotating, general-relativistic stellar model. The departures from equilibrium are analyzed into tensor spherical 
harmonics and then into complex normal modes with various mixtures of incoming and outgoing gravitational waves. 
A discussion is given of the expansion of real, physical pulsations with purely outgoing gravitational radiation in 
terms of the complex normal modes. Criteria are developed for stability against non-radial pulsations; and methods 
are devised for computing numerically the pulsation frequencies, eigenfunction, and gravitational-radiation damping 
times of the stable, real quasi-normal modes of pulsation. 



Gravitational radiation damping — Kip S. Thornc, Phys. Rev. Lett. 21, 320 (1968) 

Abstract: An analysis in presented of the emission of gravitational waves by a star in non-radial pulsation. This 
analysis, rigorous to first order in the pulsation amplitude and to all orders in G and c, exhibits damping of the 
pulsations as the waves are emitted (gravitational radiation reaction). 

Non-radial pulsation of general-relativistic stellar models. II. Properties of the grav- 
itational waves — Richard Price, Kip S. Thornc, Astrophys. J. 155, 163 (1969) 

Abstract: An analysis in presented of the properties of the gravitational waves emitted by a fully relativistic, 
non-rotating stellar model, which is pulsating in a non-radial manner. Attention is focused on quasi-normal-mode 
pulsations in which the star is excited at some initial moment of time, it radiates gravitational waves thereafter, and 
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these waves gradually damp its nearly periodic pulsations. For such quasi-normal modes a discussion is presented of 
the polarization of the gravitational waves and of the energy, momentum, and angular momentum which the waves 
transport. 

Non-radial pulsation of general-relativistic stellar models. III. Analytic and numerical 
results for neutron stars — Kip S. Thorne, Astrophys. J. 158, 1 (1969) 

Abstract: Numerical techniques axe developed for calculating the properties of the non-radial, quasi-normal pulsa- 
tions of fully relativistic stellar models. These techniques are then applied to calculate, for four realistic neutron-star 
models, the following characteristics of the lowest few quadrupole p-type modes: the pulsation periods (~ 10~^ sec), 
the gravitational-radiation damping times (~ Isec), the pulsation energies (~ 10^^ ergs), the power emitted in grav- 
itational waves (~ 10^^ ergs sec~^), and the shapes of the fluid motions. Analytic techniques are used to show that 
many, and perhaps all, of the g-type modes of a neutron star have zero frequency, if thermal pressures are negligible 
compared with nuclear pressures. 

Non-radial pulsation of general-relativistic stellar models. IV. The wecik-field limit — 

Kip S. Thorne, Astrophys. J. 158, 997 (1969) 

Abstract: In previous papers of this series techniques were developed for calculating the properties of the non-radial, 
quasi-normal pulsations of fully relativistic stellar models. In this paper those techniques are reexamined in the 
weak-field, nearly Newtonian limit, and are found to predict an emission of gravitational radiationat a rate which 
agrees with the standard weak-field approximation to general relativity. The present analysis also exhibits explicitly 
the radiation-reaction forces inside the star and shows that they sap pulsation energy from the star at the same 
rate as the waves carry it off. Following a train of ideas due to William L. Burke, the analysis is generalized to 
any weak-field, slow-motion system; and the relationship of the analysis to the post-Newtonian expansion schemes 
of Chandrasekliar and Fock is discussed. It is shown that the dominant radiation-reaction forces associated with 
each changing multipole moment can be incorporated into Newtonian theory by a simple change in the bounday 
condition for the Newtonian potential at radial infinity. 



Non-radial pulsation of general-relativistic stellar models. VI. Corrections — James R. 
Ipser, Kip S. Thorne, Astrophys. J. 181, 181 (1973) 

Abstract: The equations governing the non-radial pulsations of general-relativistic stellar modelsare shown to form 
a fourth-order system, rather than a fifth-order system as was believed previously. Incorrect statements about the 
application of bounday conditions at the stellar surface are also rectified. The errors corrected here do not invalidate 
any of the final results obtained in previous papers in this series. 

A variational principle and a stability criterion for the non-radial modes of pulsa- 
tions of stellar models in general relativity — Steven L. Detweiler, James R. Ipser, Astro- 
phys. J. 185, 685 (1973) 

Abstract: Small non-radial perturbations of spherical stellar models in general relativity are studied. A variational 
principle for the normal modes of non-radial pulsationswith / > 2 is derived. In this principle, an expression for 
the squared frequency is stationary with respect to arbitrary variations of three independent dynamical variables 
if and only if the variations are carried out about the eigenfunctions of a normal mode. The stationary value of 
the frequency is the associated eigenfrequency. For normal modes with purely outgoing radiation, the variational 
principle is modified to yield an approximate scheme for computing first the real part and then the imaginary part 
of an outgoing eigenfrequency with small imaginary part. Associated with the variational principle is an energy-like 
quantity that obeys a conservation law. It is shown that a model has no neutral, or zero-frequency, modes if the 
density decreases outwards monotonically and if the relativistic Schwarzschild discriminant is nonzero. If the dis- 
criminant is actually positive, the outgoing normal modes are stable. 



Gravitational perturbations of spherically symmetric systems. I. the exterior problem 

— Vincent Moncricf, Annals of Phys. 88, 323 (1974) 

Abstract: Gravitational perturbations of the Schwarzschild metric are treated from a point of view which is adapted, 
in a natural way, to the gauge group of the perturbed Einstein equations. The metric perturbations are explicitly 
decomposed into their gauge invariant, gauge dependent and constrained parts and a variational principle for the 
perturbation equations is derived. The Regge- Wheeler and Zerilli equations are rederived and shown to have a 
gauge invariant significance. The Hamiltonian for the perturbations is constructed and used to discuss the stability 
properties of the Schwarzschild black hole. 
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Gravitational perturbations of spherically symmetric systems. II. Perfect fluid inte- 
riors — Vincent Moncrief, Annals of Phys. 88, 343 (1974) 

Abstract: Methods developed in a previous paper on perturbations of the Schwarzschild metric are here extended 
to the treatment of perturbations of perfect fluid steUar models. The perturbations of a perfect fluid sphere are 
explicitly decomposed into their gauge invariant and gauge dependent parts and a variational principle for the per- 
turbation equations is derived. The Hamiltonian for the perturbations is constructed and a sufficient condition for 
stability against non-radial, radiative perturbations is derived from it. The stability criterion is applied to two in- 
teresting classes of stellar models, polytropic white dwarfs models and high-density neutron star cores with pressure 
proportional to energy density. 

A variational calculation of the fundamental frequencies of quadrupole pulsations of 
fluid spheres in general relativity — Steven L. Dctweiler, Astrophys. J. 197, 203 (1975) 

Abstract: A variational principle is used to calculate, in a relatively simple manner, the frequencies of quadrupole os- 
cillations of relativistic neutron star models for a variety of equations of state and central densities. Typical periods 
and damping times are on the order of milliseconds and seconds, respectively, in agreement with earlier estimates 
obtained by Thorne in a different manner. For spheres of incompressible fluids, the analysis indicates that the 
fundamental quadrupole mode becomes unstable as 2M/R attains its limiting value of 8/9. 



Torsional oscillations of neutron stars — Bonny L. Schumakcr, Kip S. Thorne, Mon. Not. 
R. Astron. Soc. 203, 457 (1983) 

Abstract: Motivated by the possibility that torsional oscillations of neutron stars may be observable in the timing 
of pulsar subpulses and/or in future gravitational-wfave detectors, this paper develops the detailed mathematical 
theory of such torsional oscillations and of the gravitational waves they emit. The oscillations are analysed using the 
formulation of first-order perturbations of a fully general relativistic spherical stellar model. All sources of damping 
are ignored except gravitational radiation reaction. The perturbations are resolved into spherical harmonics, which 
decouple from each other. For each harmonic this paper presents equations of motion, an action principle, an energy 
conservation law and a Liapunow-type proof that the oscillations are always stable. Each harmonic is then resolved 
into normal modes with outgoing gravitational waves (time dependence e*'^* with u) complex) and an eigenvalue 
problem as posed for the eigenfunctions and the eigenfrequencies u). Five methods of solving the eigenvalue problem 
axe presented; three methods axe valid in general (the method of resonances, the variational method and the method 
of energy conservation) ; one is valid on the slow- motion approximation (wavelength of waves large compared to star) 
and one is valid in the weak-gravity approximation. For stellar models with weak gravity and with radially constant 
density and shear modulus the eigenvalue problem is solved analytically. An appendix develops a general theory of 
action principles for systems with radiative boundary conditions — a theory which is then used to derive the action 
principles in the body of the paper and which could be useful for a variety of other problems involving physical 
systems coupled to radiation. 



The quadrupole oscillations of neutron stars — Lee Lindblom and Steven L. Detweiler, 
Astrophys. J. Suppl. 53, 73 (1983) 

Abstract: We investigate here the influence of the equation of state of matter (in the nuclear density regime) on the 
dynamical behavior of neutron stars. The properties of the quadrupole (/ mode) oscillations of neutron stars con- 
structed from 13 equations of state are presented. An efficient algorithm for computing the complex eigenfrequencies 
of the non-radial (p and / mode) oscillations of neutron stars is presented. The results of our computations are 
compared to relevant astrophysical observations. 



On the nonradial pulsation of general relativistic stellar models Steven L. Detweiler 
and Lee Lindblom, Astrophys. J. 292, 12 (1985) 

Abstract: We present a nonsingular fourth-order system of equations which describe the nonradial pulsations of gen- 
eral relativity stellar models. These equations represent an improvement over previous discussions in the literature 
which presented flfth-order systems or singular fourth-order systems to describe these pulsations. We also present 
the nonsingular power-series solutions to these equationsnear r = 0. 



Nonradial (/-mode oscillations of warm neutron stars — P.N. McDcrmott, H.M. Van Horn, 
J.F. SchoU, Astrophys. J. 268, 837 (1983) 

Abstract: Wc have computed adiabatic I = 1 and 2, p-, /-, and gr-mode fluid oscillations of finite-temperature neutron 
star models. We use the "relativistic Cowling approximation," defined in this paper, in which all perturbations of 
the gravitational potentials are neglected. This approximation is expected to be adequate for the g-modes and even 
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gives reasonable results for the periods of p- and /-modes. [...] For the fluid g-modc oscillations wc have studied, 
however, gravitational radiation damping is exceedingly weak, typical values of r greatly exceeding the age of the 
universe. These very long gravitational radiation damping times are a consequence of the much longer g-mode 
periods, which greatly reduce the relative power radiated, and of concentration of the oscillatory motions to the 
outermost several meters of the neutron star, which reduces the effective quadrupole moment of the oscillations. 
Damping of neutron star (/-mode oscillations will thus be dominated by other processes. The confinement of the 
p-mode eigenfunctions to the crust suggests that the effects of the nonzero shear modulus (which we have not in- 
cluded) may be quite important. 



The non-radial oscillation spectra of neutron stcirs — P.N. McDermott, C.J. Hansen, H.M. 
Van Horn, R. Buland, Astrophys. J. 297, L37 (1985) 

Abstract: We have computed non-radial oscillation spectra for a number of realistic, finite-temperature, nonrotat- 
ing, nonmagnetic neutron stars containing solid crusts. The spectra consist of p-modes (with periods If < 0.2 ms), 
g-modes trapped in the fluid "ocean" layer (IT > 100 ms), g-modes trapped in the fluid core (If > 10 s), shear (s-) 
and toroidal (t-) modes conflned to the solid crust (If < 1 ms), and interfacial (j-)modes associated with the stellar 
surface (/-mode: 11 « 0.2 — 0.4 ms), the "ocean" /crust interface (11 !5i 10 — 100 ms), and the crust/core interface 
(n » 10 - 100 ms). 



The accuracy of the quadrupole approximation for the gravitational radiation from 
pulsating stars — Eugene Balbinski, Steven Detweiler, Lee Lindblom, Bernard F. Schutz, Men. 
Not. R. Astron. Soc. 213, 553 (1985) 

Abstract: The methods for computing the damping times of the non-radial oscillations of stars due to gravita- 
tional radiation emission are compared: (i) solving the complete linearized equations of general relativity and (ii) 
using an approximate method based on the quadrupole formula. It is shown that the results of the fully relativis- 
tic calculation approach the results based on the quadrupole for sufficiently non-relativistic stars (i.e. GM/c^R 
sufficiently small) . the accuracy of the approximation method is investigated by considering polytropic stars with a 
range of polytropic indices. The approximation is worse for larger polytropic indices, i.e. for greater condensation. 
Wc conclude the when Newtonian models are used to approximate relativistic ones, models with the same value of 
GM/c^R (rather than, say, central density) should be compared. 

Nonradial oscillations of neutron stars — P.N. McDermott, H.M. Van Horn, C.J. Hansen, 
Astrophys. J. 325, 725 (1988) 

Abstract: Finite-temperature neutron star models with a fluid core, solid crust, and thin surface fluid "ocean" 
are capable of sustaining a broad diversity of normal modes of oscillations. We have performed linear, adiabatic, 
Newtonian, non-radial pulsation analyses of such stars, including the clastic effects of the neutron star crust. [...] 
Damping mechanisms we have investigated include gravitational radiation damping, neutrino emission damping, 
electromagnetic radiation from an oscillating stellar magnetic fleld, nonadiabatic effects, and internal friction and 
viscosity. Which of these mechanisms dominates depends upon the type of mode, the spherical harmonic index I, 
the equilibrium model, and other parameters. Finally, we suggest that surface g-modes may be excited during the 
thermonuclear outbursts in X-ray bursters and may lead to observable quasi-periodic oscillations. 



Normal modes of a model radiating system — Kostas. D. Kokkotas and Bernard F. Schutz, 

Gen. Relativity Gravitation 18, 913 (1986) 

Abstract: In order to gain insight into normal modes of realistic radiating systems, we study the simple model 
problem of a finite string and a semi-inflnite string coupled by a spring. As expected there is a family of modes 
which are basically the modes of the finite string slowly damped by the "radiation" of energy to infinity on the 
semi-infinite string. But we also study another family of modes, found by Dyson in a different model problem, which 
are strongly damped modes of the semi-inflnite string itself. These may be analogous to the modes of black holes, 
and they are likely to be present in relativistic stars as well. The question of whether the instability in these modes 
which Dyson found is present in realistic stars remains open. 



Stellar resonant oscillations coupled to gravitational waves — Yasufumi Kojima, Prog. 
Theo. Phys. 77, 297 (1987) 

Abstract: Gravitational waves emitted by a particle moving in a circular orbit around a stellar model are examined. 
First-order perturbation equations are solved both inside and outside the star. A resonant oscillation of the star oc- 
curs when the frequency of the gravitational waves produced by the particle coincides with the quasi-normal mode of 
the star. The amplitude of the gravitational wave at the resonance exhibits a sharp intrinsic peak depending on the 
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equation of state used in the equilibrium configuration. 



Two families of normal modes in relativistic stars — Yasufumi Kojima, Prog. Theo. Phys. 
79, 665 (1988) 

Abstract: Normal modes of relativistic stars are examined in order to study the systems consisting of a star and 
gravitational waves. It is found that there are two classes of normal modes in such a radiating system. They are 
extremely different in the decay rates of the complex eigenfrcquencies. One class is a well-known family of weakly 
damped modes and the other is a family of strongly damped modes, which are discovered for the first time in 
numerical calculations for the relativistic stellar oscillations. The weakly damped modes represent the slow decay 
of the stellar oscillations due to gravitational radiation. As the system becomes more relativistic, the emission 
rate of gravitational waves increases. As a result the decay rates of these modes increase. On the other hand, the 
strongly damped modes essentially represent rapid decay of the gravitational waves penetrated into the star because 
displacements of matter inside the star are found to be very small. The decay rates of these modes decrease as the 
system becomes more relativistic. 



The effect of viscosity on neutron star oscillations — Curt Cutler, Lee Lindblom, Astrophys. 
J. 314, 234 (1987) 

Abstract: In this paper we compute the rate at which neutron star oscillations arc damped by the presence of 
viscosity and thermal conductivity in the neutron matter. In our computation we use fully relativistic equations to 
describe the neutron star oscillations, we use the best available expressions for the dissipation coefficients in neutron 
matter, and we approximate the effects of superfluidity on the dissipation mechanisms. We compute damping times 
for a range of neutron star masses and temperatures, and we use a variety of equations of state to describe the 
highest density portion of the neutron matter.The time scales computed here are used to improve previous estimates 
of the maximum rotation rate of neutron stars. 



Non-radial pulsations of neutron stars with a crust — Lee Samuel Finn, Mon. Not. R. 
Astron. Soc. 245, 82 (1990) 

Abstract: A fully relativistic treatment of non-radial pulsational modes of spherical stars with solid crusts is pre- 
sented. This work complements Schumaker & Thome's relativistic treatment of torsional modes. The stain tensor 
for arbitrary even-parity elastic deformations of a spherical star is derived. A Hookean law is assumed to relate the 
stress to the strain, and the pulsational equations and their bounday conditions are derived from the Einstein field 
equations and the equations of motion. The Einstein field equations and equations of motion arc tlicu used to find 
the Lagrangian density and an action principle for the pulsational problem. This action principle may prove useful 
in future studies of the stability of non-radial pulsational modes for stars with solid crust. 



On the non-radial oscillations of a star — Subrahmanyan Chandrasekhar, Valeria Ferrari, 
Proc. R. Soc. London, Ser. A 432, 247 (1991) 

Abstract: A complete theory of the non-radial oscillations of a static spherically symmetric distribution of matter, 
described in terms of an energy density and an isotropic pressure, is developed, ab initio, on the premise that the 
oscillations are excited by incident gravitational waves. The equations, as formulated, enable the decoupling of the 
equations governing the perturbations of the metric of the space-time from the equations governing the hydrody- 
namical variables. This decoupling of the equations reduces the problem of determining the complex characteristic 
frequencies of the quasi-normal modes of the non-radial oscillations to a problem in the scattering of incident grav- 
itational waves by the curvature of the space-time and the matter content of the source acting as a potential. The 
present paper is restricted [..] to the case when the underlying equation of state is barotropic. The algorism de- 
veloped for the determination of the quasi-normal modes is directly confirmed by comparison with an independent 
evaluation by the extant alternative algorism. Both polar and axial perturbations are considered. Dipole oscilla- 
tions [...] are also treated as a particular simple special case. [...] 

On the non-radial oscillations of a star. II Further amplifications — Subrahmanyan 
Chandrasekhar, Valeria Ferrari, Roland Winston Proc. R. Soc. London, Ser. A 434, 635 (1991) 

Abstract: The two algorisms that one uses for determining the complex characteristic frequencies, cr = cro-hicrj, 
belonging to the quasi-normal modes of oscillations of a star, are shown to be equivalent. In the first algorism, one 
searches directly in the complex cr-plane to satisfy the requirements that at infinity there are only outgoing waves 
(and no iiiconiiiig waves), in the second algorism, one restricts oneself to real crs and selects the solution with the 
minimum fiux of gravitational radiation at infinity. 
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On the non-radial oscillations of a star. Ill A reconsideration of the axial modes — 
Subrahmanyan Chandrasekhar, Valeria Ferrari, Proc. R. Soc. London, Ser. A 434, 449 (1991) 
Abstract: It is shown that for stars with radii in the range 2.1hGM j (? < R < ca,. dGM/c^, quasi-normal axial 
modes of oscillations are possible. These modes are explicitly evaluated for stellar models of uniform density. 

On the non-radial oscillations of a star. IV An application of the theory of Regge 
poles — Subrahmanyan Chandrasekhar, Valeria Ferrari, Proc. R. Soc. London, Ser. A 437, 133 
(1992) 

Abstract: It is shown how the flow of energy in the form of gravitational radiation, through a star in non-radial os- 
cillations, can be determined by a suitable adaption of Rcgge's theory of potential scattering in the quantum theory. 
While the resonant scattering of axial gravitational waves can be treated by Regge's theory, essentially, in its original 
form, the treatment of resonant scattering of polar gravitational waves requires a generalization of the theory beyond 
its standard range and the incorporation of a 'flux integral' derived from independent considerations. Illustrative 
examples of the applications of the theory are provided. 



Nonradial pulsations of stellar models in general relativity — James R. Ipser, Richard H. 
Price Phys. Rev. D 43, 1768 (1991) 

Abstract: In the Newtonian theory, a method was developed recently that enables one to solve the previously in- 
tractable problem of obtaining solutions to the eigenequations for the normal modes of rotating stellar models. In 
this paper, a closely related method is used to treat the even-parity non-radial modes of spherical stellar mod- 
els in general relativity. In this treatment of generally nonbarotropic perturbations, the Lagrangian displacement is 
eliminated by solving for it in terms of a certain fluid perturbation quantity 5h (the enthalpy perturbation in the 
barotropic case) and the Regge- Wheeler metric perturbations K and Hq. One is then led to an explicetly fourth- 
order system, consisting of two secoud-oidci equations for either pair K and 5h or K and Ho- The fourth-order 
system for K and Hq, from which all fluid perturbation variables have been eliminated, is relatively simple and 
provides a new way to compute quasi-normal non-radial modes and to analyze the interaction of a stellar model with 
incident gravitational wave. The analysis provide insight into the connection between the Newtonian and relativis- 
tic theories of stellar pulsations. The relationship of the analyses to the recent work of Chandrasekhar and Ferrari 
is discussed. 



Relation of gauge formalisms for pulsations of general-relativistic stellar models — 

Richard H. Price, James R. Ipser, Phys. Rev. D 44, 307 (1991) 

Abstract: 

Initial-boundary value problem for the spherically symmetric Einstein equations for 
a perfect fluid Saskia Kind, Jiirgen Ehlers, Class. Quantum Grav. 10, 2123 (1993) 

Abstract: It is shown that for a given spherically symmetric distribution of a perfect fluid on a spaeelike hypersur- 
face with boundary and a given, time-dependent boundary pressure, there exists a unique, local-in-time solution of 
the Einstein equations. A Schwarzschild spacetime can be attached to the fluid body if and only if the boundary 
pressure vanishes. We assume a smooth equation of state for which the density and the speed of sound remain 
positive for vanishing pressure. 



Relativistic stellar oscillations treated as an initial value problem — Saskia Kind, Jiirgen 
Ehlers, Bernd G. Schmidt, Class. Quantum Grav. 10, 2137 (1993) 

Abstract: The linearized Einstein equations for a static, spherically symmetric fluid ball and its empty surroundings 
are considered. It is shown that, given initial data obeying the constraints, there exists a unique solution, which 
describes the motion of the perturbed fluid and the gravitational waves propagating inside and outside the fluid 
ball. The physical junction conditions for the boundary of the ball suffice to determine the evolution inside and 
outside of the ball in terms of initial values. The equation of state is assumed smooth and such that the density 
and the speed of sound remain positive for vanishing pressure. 



VK-modes: a new family of normal modes of pulsating relativistic stars — Kostas D. 
Kokkotas, Bernard F. Schutz, Mon. Not. R. Astron. Soc. 255, 119 (1992) 

Abstract: We demonstrate explicitly the existence of a new family of outgoing-wave normal modes of pulsating rel- 
ativistic stars, the first such family known that has no analogue in Newtonian stars. These modes were discovered 
by the authors in a toy model, where they were called strongly damped normal modes. Kojima then found the 
first examples of these modes in realistic spherical polytropic stellar models. Here we give a number of arguments 
that demonstrate the existence of this family unequivocally, and we calculate a large number of eigenfrequencies. 
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Physically, the modes arise directly from the coupling of the fluid oscillations of the star to the gravitational wave os- 
cillationsof the space-time metric. Previously studied modes of relativistic stars have been close analogues of modes 
of Newtonian stars, where the eoupling to gravitational waves mainly generated a small imaginary part to the fre- 
quency. Such modes can be classified using the Newtonian classes: /-, p-, g- and r-modes. The present modes, by 
contrast, have no Newtonian analogue. They are primarily oscillations of the space-time metric in which the fluid 
hardly vibrates at all. We christen them w-modes (gravitational-wave modes). These modes are strongly damped, 
being characterized by complex frequencies with unusually large imaginary parts, comparable to their real parts. 
We calculate a sequence of Z = 2 modes for a number of spherical polytropic stellar models. An interesting feature 
of to-modes is that the lowest order mode of each sequence has a frequency similar to that of lowest order mode of a 
black hole. For higher modes, the spectrum diverges from the black hole spectrum, but shows remarkable similarity 
to that of the strongly damped modes of the toy problem. As carriers of gravitational wave information, w-modes 
may be important and observable in the burst of gravitational radiationthat follows the formation of a neutron star. 
They should also be essential in solving the problem of the completeness of the outgoing-wave normal modes of 
radiating systems. 



Nonradial oscillations of neutron stars: A new branch of strongly damped normal 
modes — Markus Leins, Hans-Peter NoUert, Michael H. Soffel, Phys. Rev. D 48, 3467 (1993) 

Abstract: We study small, non-radial oscillationsof neutron stars in a general relativistic perturbation treatment, 
considering different values for the central density of the star. We adapt two techniques used previously for the 
determination of quasi-normal mode of black holes, allowing us for the first time to determine without approxi- 
mation solutions which satisfy purely outgoing bounday conditions at spatial infinity. We confirm the existence of 
strongly damped complex normal modes found by Kokkotas and Schutz (w-modes). In addition, we identify a new 
branch of strongly damped modes (w //-modes) . Our new modes are much more similar to quasi- normal modes of 
black holes than the w-modes known before. 



Axial modes for relativistic stars — Kostas D. Kokkotas, Men. Not. R. Astron. Soc. 268, 
1015 (1994) + Erratum, Men. Not. R. Astron. Soc. 277, 1599 (1995) 

Abstract: For the class of stellar oscillations named ocZcZ-parity of axial modes, Chandrasekliar & Ferrari have found 
that, under appropriate conditions, there exists a spectrum with frequencies describing quasi-normal oscillations hav- 
ing long damping times. In this work, we use different (and more appropriate for this kinf od problem) numerical 
methods to show the existence of an additional part of the axial mode spectrum with quasi-normal mode frequen- 
cies having short damping times, for precisely the same stellar models that Chandrasekhar &; Ferrari have considered. 
This part of the axial modes spectrum is similar to the w-modes found earlier by Kokkotas & Schutz. 



A new numerical approach to the oscillation modes of relativistic stars Nils Andcrsson, 
Kostas D. Kokkotas, Bernard F. Schutz, Mon. Not. R. Astron. Soc. 274, 1039 (1995) 

Abstract: The oscillation modes of a simple polytropic stellar model are studied. Using a new rmmcrieal approach 
(based on integration for complex coordinates) to the problem for the stellar exterior we have computed the eigen- 
frequencies of the highly damped w-modes. The results obtained agree well with the recent ones of Leins, NoUert & 
Soffel. Specifically, we are able to explain why several modes in this regime of the complex frequency plane could not 
be identified within the WBK approach of Kokkotas & Schutz. Furthermore, we have established that the 'kink' that 
was a prominent feature of the spectra of Kokkotas & Schutz, but which did not appear in the results of Leins et al., 
was a numerical artefact. Using our new numerical code we are also able to compute, for the first time, several of the 
slowly damped (p)-modes for the considered stellar models. For very compact stars, we find, somewhat surprisingly, 
that the damping of these modes does not decrease monotonically as one proceeds to higher oscillation frequencies. 
The existence of low-order modes that damp away much faster that anticipated may have implications for questions 
regarding stellar stability and the lifetime of gravitational-wave sources. The present results illustrate the accuracy 
and reliability of the complex-coordinate method and indicate that the method could prove to be of great use also 
in problems involving rotating stars. There is no apparent reason why the complex-coordinate approach should not 
extend to rotating stars, whereas it is accepted that all previous methods will fail to do so. 



Two simple models for gravitational wave modes of compact stars — Nils Andersson 

Gen. Relativity Gravitation 28, 1433 (1996) 

Abstract: Two simple model problems relevant for the gravitational wave modes of relativistic stars are discussed. It 
is shown that the entire mode-spectrum can be obtained if one considers the modes as arising because of the trapping 
of gravitational waves by the spaeetime curvature. The stellar fiuid need play no dynamical role. Furthermore, it is 
shown that two distinct families of gravitational wave modes exist. The first corresponds to waves trapped inside 
the star, while the second is similar to acoustic waves scattered off a hard sphere. An infinite number of the first 
kind of modes exist, but the latter family will only have a few members. 
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On the oscillation spectra of ultracompact stars: An extensive survey of gravitational- 
wave modes — Nils Andersson, Yasufumi Kojima, Kostas D. Kokkotas Astrophys. J. 462, 855 

(1996) 

Abstract: An extensive survey of gravitational- wave modes for uiiiforiii density stars is presented. The study covers 
ranging in compactness from R/M = 100 to the Umit of stabihty in general relativity: R/M = 9/4. We establish 
that polar and axial gravitational-wave modes exist for all these stellar models. Moreover, there are two distinct 
families of both axial and polar modes. We discuss the physics of these modes and argue that one family is primarily 
associated with the interior of the star, while the second family is mainly associated with the stellar surface. We also 
show that the problem of axial perturbations has all the essential features of the polar problem as far as gravitational 
waves are concerned. This means that the axial problem is much more important that has previously been assumed. 
We also find some surprising features, auch as avoided crossings between the polar gravitational-wave modes and 
the Kelvin /-mode as the star becomes very compact. This seems to suggest that the /-mode should be considered 
on equal footing with the polar to-modes for ultracompact stars. All modes may have the main character of trapped 
modes inside the curvature potential barrier for R/M < 3. 



Space-time modes of relativistic stars — Nils Andersson, Kostas D. Kokkotas, Bernard F. 

Schutz, Mon. Not. R. Astron. Soc. 280, 1230 (1996) 

Abstract: The problem of relativistic stellar pulsations is studied in a somewhat ad hoc approximation that ignores 
all fluid motions. This Inverse Cowling Approximation (ICA) is motivated by two observations. (1) For highly 
damped (ly-mode) oscillations the fluid plays very little role. (2) If the fluid motion is neglected, the problem for 
polar oscillation modes becomes similar to that for axial modes. Using the ICA, we find a polar mode spectrum that 
has all features of the -iii-mode spectrum for the full problem. Moreover, in the limit of superdense stars, we find 
the ICA spectrum to be qualitatively similar to that of the axial modes. These results clearly show the importance 
of general relativity for the pulsation modes of compact stars, and that there are modes whose existence does not 
depend on motions of the fluid at all, namely pure 'space-time' modes. 



Gravitational waves and pulsating stcirs: What can we learn from future observations? 

— Nils Andersson, Kostas D. Kokkotas, Phys. Rev. Lett. 77, 4134 (1996) 

Abstract: We present new results for pulsating stars in general relativity. First, we show that the so-called gravita- 
tional modes of a neutron star can be excited when a gravitational wave impinges on the star. Numerical simulations 
suggest that the modes may be astrophysically relevant, and we discuss whether they will be observable with future 
gravitational wave detectors. We also discuss how such observations could lead to estimates of both the radius and 
the mass of a neutron star, and thus put constraints on the nuclear equation of state. 



Relativistic stellar pulsations with near-zone bounday conditions — Lee Lindblom, Gre- 
gory Mendell and James R. Ipser, Phys. Rev. D 56, 2118 (1997) 

Abstract: A new method is presented here for evaluating approximately the pulsation modes of relativistic stellar 
models. This approximation relies on the fact that gravitational radiation influences these modes only on time scales 
that are much longer that the basic hydrodyiiamical time scale of the system. This makes it possible to imjjose the 
bounday conditions on the gravitational potentials at the surface of the star rather than in the asymptotic wace 
zone of the gravitational field. This approximation is tested here by predicting the frequencies of the outgoing non- 
radial hydrodynamie modes of nonrotating stars. The real parts of the frequencies are determined with an accuracy 
that is better than our knowledge if the exact frequencies (about 0.01%) except in the most relativistic models 
where it decreases to about 0.1%. The imaginary parts of the frequencies are determined with an accuracy of 
approximately M/ R, where M is the mass and R is the radius of the star in question. 

Gravitational waves from pulsating stars: Evolving the perturbation equations for a 
relativistic star Gabriclk; Allen, Nils Andersson, Kostas D. Kokkotas, Bernard, F. Schutz 
Phys. Rev. D 58, 124012 (1998) 

Abstract: We consider the perturbations of a relativistic star as an initial value problem. Having discussed the 
formulation of the problem (the perturbation equations and the appropriate bounday conditions at the center and 
the surface of the star) in detail we evolve the equations numerically from several different sets of initial data. In 
all considered cases we find that the resulting gravitational waves carry the signature of several of the star's pulsa- 
tion modes. Typically, the fluid /-mode, the flrst two p-modes and the slowest damped gravitational ui-mode are 
present in the signal. This indicates that the pulsation modes may be an interesting source for detectable gravita- 
tional waves from colliding neutron stars or supernovae. We also survey the literature and find several indications of 
mode presence in numerical simulations of rotating core collapse and coalescing neutron stars. If such mode-signals 
can be detected by future gravitational wave antennae one can hope to infer detailed information about neutron 
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stars. Since a perturbation evolution should adequately describe the late time behavior of a dynamically excited 
neutron star, the present work can also be used as a bench-mark test for future nonlinear simulations. 



Pulsation modes for increasingly relativistic polytropes — Nils Andersson, Kostas D. 
Kokkotas, Mon. Not. R. Astron. Soc. 297, 492 (1998) 

Abstract: We present the results of a numerical study of the fluid /-, p- and the gravitational lo-modes for increas- 
ingly relativistic non-radial polytropes. The results for /- and lo-modes are in good agreement with previous data 
for uniform density stars, which supports an understanding of the nature of the gravitational wave modes based on 
the uniform density data. We show that the p-modes can become extremely long-lived for some relativistic stars. 
This effect is attributed to the change in the perturbed density as the star becomes more compact. 



Towards gravitational wave asteroseismology — Nils Andersson, Kostas D. Kokkotas, Mon. 
Not. R. Astron. Soc. 299, 1059 (1998) 

Abstract: We present new results for pulsating neutron stars, we have calculated the eigenfrequencies of the modes 
that one would expect to be the most important gravitational wave sources; the fundamental fluid /-mode, the 
first pressure p-mode and the first gravitational wave ui-mode, for twelve realistic equations of state. From these 
numerical data we have inferred a set of 'empirical relations' between the mode frequencies and the parameters of 
the star (the radius R and the mass M). Some of these relations prove to be surprisingly rebust, and we show 
how they can be used to extract the details of the star from observed modes. The results indicate that, should the 
various pulsation modes be detected by the new generation of gravitational wave detectors that come online in a 
few years, the mass and the radius of neutron stars can be deduced with errors no larger that a few per cent. 



The imprint of the equation of state on the axial w-modes o f oscillating neutron stars 

— Omar Benhar, Emanuclc Berti, Valeria Ferrari, gr-qc/9901037 

Abstract: We discuss the dependence of the pulsation frequencies of the axial quasi-normal modes of a non-radial neu- 
tron star upon the equation of state describing the star interior. The continued fraction method has been used to 
compute the complex frequencies for a set of equations of state based on different physical assumptions and spanning 
a wide range of stiffness. The numerical results show that the detection of axial gravitational waves would allow 
to discriminate between models underlying the different equations of state, thus providing relevant information on 
both the structure of neutron star matter and the nature of the hadronic interactions. 



The inverse problem for pulsating neutron stars: A "fingerprint analysis " for the 
supranuclear e quation of state — Kostas D. Kokkotas, T.A. Apostolatos and Nils Andersson, 
gr-qc/9901072| 

Abstract: We study the problem of detecting, and infering astrophysical information from, gravitational waves from 
a pulsating neutron star. We show that the fluid / and p-modes, as well as the gravitational wave ui-modes may be 
detectable from sources in our own galaxy, and investigate how accurately the frequencies and damping rates of these 
modes can be infered from a noisy gravitational wave data stream. Based on the conclusions of this discussion we 
propose a strategy for revealing the supranuclear equation of state using the neutron star fingerprints: the observed 
frequencies of an / and a p-mode. We also discuss how well the source can be located in the sky using observations 
with several detectors. 



Construction of astrophysical initial data for perturbations of relativistic stars — Nils 
Andersson. Kost as D. Kokkotas, Pablo Laguna, Phillippos Papadopoulos and Michael S. Sipior, 
gr-qc/9904059| 

Abstract: We develop a framework for constructing initial data sets for perturbations about spherically symmetric 
matter distributions. This framework faciliates setting initial data representing astrophysical sources of gravita- 
tional radiation involving relativistic stars. The procedure is based on Lichnerowicz- York's conformal approach to 
solve the constraints in Einstein's equations. The correspondence of these initial data sets in terms of the standard 
gauge perturbation variables in the Regge- Wheeler perturbation variables is established, and examples of initial 
data sets of merging neutron starsunder the close-limit approximation are pesented. 
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Papers on Particles orbiting Black Holes and Neutron Stars 



Gravitational field of a particle falling in Schwarzschildgeometry analyzed in tensor 
harmonics Frank J. Zcrilli, Phys. Rev. D 2, 2141 (1970) 

Abstract: We arc concerned with the pulse of gravitational radiation given off when a star falls into a "blacl; hole" 
near the center of out galaxy. We look at the problem of a small particle falling in a Schwarzschild background 
("black hole") and examine its spectrum in the high-frequency limit. In formulating the problem it is essential to 
pose correct bounday conditions: gravitational radiation not only escaping to infinity but also dissappearing down 
the hole. We have examined the problem in the approximation of linear perturbations from a Schwarzschild back- 
ground geometry, utilizing the decomposition into tensor spherical harmonics given by Regge and Wheeler (1957) 
and by Mathews (1962). The falling particle contributes a 5-function source term (geodesic motion in the background 
Schwarzschild geometry) which is also decomposed into tensor harmonics, each of which "drives" the corresponding 
perturbation harmonic. The power spectrum radiated in infinity is given in the high frequency approximation in 
terms of the traceless transverse tensor harmonics called "electric" and "magnetic" by Mathews. 



Stellar resonant oscillations coupled to gravitational waves — Yasufumi Kojima, Prog. 
Theo. Phys. 77, 297 (1987) 

Abstract: Gravitational waves emitted by a particle moving in a circular orbit around a stellar model are examined. 
First-order perturbation equations are solved both inside and outside the star. A resonant oscillation of the star 
occurs when the frequency of the gravitational waves produced by the particle coincides with the quasi-normal 
mode of the star. The amplitude of the gravitational wave at the resonance exhibits a sharp intrinsic peak depend- 
ing on the equation of state used in the equilibrium configuration. 

Gravitational radiation from a particle in circular orbit around a black hole. I. Ana- 
lytic results for the nonrotating case — Eric Poisson, Phys. Rev. D 47, 1497 (1993) 
Abstract: [...] This paper is the first in a series that will carry out detailed calculations relevant to [...] binaries, in 
the case where one body is a small-mass black hole of neutron star and the other is a much more massive black hole, 
and the orbit is circular (aside from its gradual inspiral). These papers will focus primarily on the emitted waveforms 
and especially their phasing — which is crucial for extraction of information from the detectors' measurements. This 
first paper is restricted to the case where the massive black hole is non-radial. The paper begins by bringing the al- 
ready well-developed formalisms for computing waveforms (the "Regge- Wheeler" and "Teukolsky" formalisms) into 
a combined form that is particularly well suited both for high accuracy numerical calculations (carried out in paper 
II), and for analytic calculations (this paper). [...] It is shown that propagation of the waves through the intermediate 
zone (which connects the near zone and the wave zone) distorts the waveforms and changes their power (and hence 
phasing), at post^/^-Newtonian order, in ways that have not previously been computed - except abstractly and 
nonconcretely as formal "tail terms" in the waves. It is demonstrated that these post^/^-Newtonian corrections will 
be of considerable importance for the extraction of information from the waveforms that LIGO/ VIRGO expects to 
measure. 



Gravitational radiation from a particle in circular orbit cU"ound a black hole. I. Nu- 
merical results for the nonrotating case Curt Cutler, Lee Samuel Finn, Eric Poisson, 
Gerald Jay Sussman, Phys. Rev. D 47, 1511 (1993) 

Abstract: [...] We consider the particular case of a compact object (e.g., either a neutron star or a stellar mass 
black hole) in a circular orbit about a much more massive Schwarschild black hole. In this limit, the gravita- 
tional radiation luminosity can be calculated by integrating the Teukolsky equation. Numerical integration us used 
to obtain accurate estimates of the luminosity dE/dt as a function of the orbital radius ro- These estimated are 
fitted to a post-Newtonian expansion of the form dE/dt = {dE/dt)p; ife^:*, where (dE/dt)ff is the standard 
quadrupolc formula expression and x = {M/tq)^/^ . From our fits wc obtain values for the expansion coefficients 
aj; up through order . While our results arc in excellent agreement with low-order post-Newtonian calculations, 
we find that the post-Newtonian expansion converges slowly. Corrections beyond may be needed to achieve the 
desired accuracy for the construction of template waveforms. 

Gravitational radiation emitted when a mass falls onto a compact star — A. Borelli, 
Nuovo Cimento B 112, 225 (1997) 
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Abstract: In this paper we study the energy spectrum related to the axial perturbations of a compact star when a 
particle falls spiralling onto it. We find that both slowly-damped quasi-normal modes and strongly damped ui-modes 
are excited, and that a part of the energy in the process is associated to these lu-modes. Our analysis will show 
substantial difference between the energy spectrum of compact stars and black holes. 



Stellar pulsa tions excited by a scattered mass — Valeria Ferrari, L. Gualtieri, A. Borelli, 
gr-qc/9901060| 

Abstract: We compute the energy spectra of the gravitational signals emitted when a mass mg is scattered by the 
gravitational field of a star of mass M 2> mo- We show that, unlike black holes in similar processes, the quasi-normal 
modes of the star are excited, and that the amount of energy emitted in these modes depends on how close the 
exciting mass cna get to the star. 



Gravitational waves from a test particle scattered by a neutron star: Axial mode case 

— Kazuhiro Tominaga, Motoyuki Saijo, Phys. Rev. D 60, 024004 (1999) 

Abstract: Using a metric perturbation method, we study gravitational waves from a test particle scattered by a 
spherically symmetric relativistic star, we calculate the energy spectrum and the waveform of gravitational waves for 
axial modes. Since metric perturbations in axial modes do not couple to the matter fluid of the star, emitted waves 
for a neutron star show only one peak in the spectrum, which corresponds to the orbital frequency at the turning 
point, where the gravitational field is strongest. However, for an ultracompact star (the radius R < 3M), another 
type of resonant periodic peak appears in the spectrum. This is just because of an excitation by a scattered particle 
of axial quasi-normal modes, which were found by Chandrasekhar and Ferrari. This excitation comes from the 
existence of the potential minimum inside the star. We also find for an ultracompact star many small periodic 
peaks at the frequency region beyond the maximum of the potential, which would be due to a resonance of two 
waves reflected by two potential barriers (Regge- Wheeler type and one at the center of the star). Such resonant 
peaks appear neither for a normal neutron star nor for a Schwarzschild black hole. Consequently, even if we analyze 
the energy spectrum of gravitational waves only for axial modes, it would be possible to distinguish between an 
ultracompact star and a normal neutron star (or a Schwarzschild black hole). 



Excitation of the od d-parity quas i-normal modes of compact objects — Zeferino Andrade 
and Richard H. Price, |gr-qc/990206^ 

Abstract: The gravitational radiation generated by a particle in a close unbounded orbit around a neutron star is 
computed as a means to study the importance of the w modes of the neutron star. For simplicity, attention is 
restricted to odd-parity ("axial") modes which so not couple to the neutron star's fluid modes. We find that for 
realistic neutron star models, particles in unbounded orbits only weakly excite the w modes; we conjecture that this 
is also the case for astrophysically interesting sources of neutron star perturbations. We also find that for cases in 
which there is significant excitation of quadrupole w modes, there is comparable excitation of higher multipoles. 
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